Disorder-Induced Critical Phenomena in Hysteresis: 
A Numerical Scaling Analysis 
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Experimental systems with a first order phase transition will often exhibit hysteresis when out of 
equilibrium. If defects are present, the hysteresis loop can have different shapes: with small disorder 
the hysteresis loop has a macroscopic jump, while for large disorder the hysteresis loop is smooth. 
The transition between these two shapes is critical, with diverging length scales and power laws. 
We simulate such a system with the zero temperature random field Ising model, in 2, 3, 4, 5, 7, 
and 9 dimensions, with systems of up to 1000 3 spins, and find the critical exponents from scaling 
collapses of several measurements. The numerical results agree well with the analytical predictions 
from a renormalization group calculation |l4[. 
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I. INTRODUCTION 

The increased interest in real materials in condensed 
matter physics has brought disordered systems into the 
spotlight. Dirt changes the free energy landscape of a 
system, and can introduce metastable states with large 
energy barriers 0. This can lead to extremely slow re- 
laxation towards the equilibrium state. On long length 
scales and practical time scales, a system driven by an 
external field will move from one metastable local free- 
energy minimum to the next. The equilibrium, global 
free energy minimum and the thermal fluctuations that 
drive the system toward it, are in this case irrelevant. 
The state of the system will instead depend on its his- 
tory. 

The motion from one local minima to the next is a col- 
lective process involving many local (magnetic) domains 
in a local region - an avalanche. In magnetic materi- 
als, as the external magnetic field H is changed contin- 
uously, these avalanches lead to the magnetic noise: the 
Barkhausen effect This effect can be picked up as 

voltage pulses in a coil surrounding the magnet. The 
distribution of pulse (avalanche) sizes is found |^-|(| to 
follow a power law with a cutoff after a few decades, and 
was interpreted by some Q to be an example of self- 
organized criticality [Q. (In SOC, a system organizes 
itself into a critical state without the need to tune an ex- 
ternal parameter.) Other systems can exhibit avalanches 
as well. Several examples where disorder may play a 
part are: superconducting vortex line avalanches ||, re- 
sistance avalanches in superconducting films JO], and cap- 
illary condensation of helium in Nuclepore pi . 

The history dependence of the state of the system 
leads to hysteresis. Experiments with magnetic tapes 
PI have shown that the shape of the hysteresis curve 
changes with the annealing temperature. The hystere- 
sis curve goes from smooth to discontinuous as the an- 
nealing temperature is increased. This transition can be 
explained in terms of a plain old critical point with two 



tunable parameters: the annealing temperature and the 
external field. At the critical temperature and field, the 
correlation length diverges, and the distribution of pulse 
(avalanche) sizes follows a power law. 

We have argued earlier |T3] that the Barkhausen noise 
experiments can be quantitatively explained by a model 
]i~3 | with two tunable parameters (external field and dis- 
order), which exhibits universal, non-equilibrium collec- 
tive behavior. The model is athermal and incorporates 
collective behavior through nearest neighbor interactions. 
The role of dirt or disorder, as we call it, is played by 
random fields. This paper presents the results and con- 
clusions of a large scale simulation of that model: the 
non-equilibrium zero-temperature Random Field Ising 
Model (RFIM), with a deterministic dynamics. The re- 
sults compare very well to our e expansion [ |T4]|l5| ] , and to 
experiments in Barkhausen noise fll2||. A more detailed 



comparison to experimental systems is in process 16 1. 

The paper is divided as follows. Section II quickly 
reviews the model. Section III explains the simulation 
method that we use. Section IV explains the data anal- 
ysis and shows results for the simulation in 2, 3, 4, and 
5 dimensions, as well as in mean field. Section V gives a 
comparison between the simulation and the e expansion 
exponents, and a comparison between the shape of the 
magnetization curves in 5, 7, and 9 dimensions, and the 
predicted shape from the e expansion. Section VI sum- 
marizes the results. This is followed by three appendices 
that cover derivations that were omitted in the text for 
continuity. 



II. THE MODEL 

To model the long-range, far from equilibrium, collec- 
tive behavior mentioned in the previous section, we define 
]i~3[ spins Si on a hypercubic lattice, which can take two 
values: s,; = ±1. The spins interact ferromagnetically 
with their nearest neighbors with a strength , and are 
sitting in a uniform magnetic field H (which is directed 
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along the spins). Dirt is simulated by a random field hi, 
associated with each site of the lattice, which is given by 
a gaussian distribution function p(hi): 
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of width proportional to R which we call the disorder 
parameter, or just disorder. The hamiltonian is then 



Jij Si Sj 



(2) 
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For the analytic calculation, as well as the simulation, we 
have set the interaction between the spins to be indepen- 
dent of the spins and equal to one for nearest neighbors, 
Jij = J = 1, and zero otherwise. 

The dynamics is deterministic, and is defined such that 
a spin Si will flip only when its local effective field : 



■ eff 



H + hi 



(3) 



changes sign. All the spins start pointing down (s, = — 1 
for all i). As the field is adiabatically increased, a spin 
will flip. Due to the nearest neighbor interaction, a 
flipped spin will push a neighbor to flip, which in turn 
might push another neighbor, and so on, thereby gener- 
ating an avalanche of spin flips. During each avalanche, 
the external field is kept constant. For large disorders, 
the distribution of random fields is wide, and spins will 
tend to flip independently of each other. Only small 
avalanches will exist, and the magnetization curve will 
be smooth. On the other hand, a small disorder im- 
plies a narrow random field distribution which allows 
larger avalanches to occur. As the disorder is lowered, 
at the disorder R = R c and field H = H c , an infinite 
avalanche in the thermodynamic system will occur for 
the first time, and the magnetization curve will show a 
discontinuity. Near R c and H c , we find critical scaling 
behavior and avalanches of all sizes. Therefore, the sys- 
tem has two tunable parameters: the external field H 
and the disorder R. We found from the mean field cal- 
culation jl4 15 and the simulation that a discontinuity 
in the magnetization exists for disorders R < i? c , at the 
field H C (R) > H C {R C ), but that only at (R C ,H C ), do we 
have critical behavior. For finite size systems of length L, 
the transition occurs at the disorder R e Jf [L) near which 
avalanches first begin to span the system in one of the d 
dimensions (spanning avalanches). The effective critical 
disorder R%" (L) is larger than i? c , and R e Jf [L) — > R c as 
L — > oo. 



III. ALGORITHM 

There are several methods that can be used to simu- 
late the above model. The simplest but most time and 



space (memory) consuming method starts by assigning a 
random field to each spin on the hypercubic lattice. At 
the beginning of the simulation, all the spins are point- 
ing down. The external field H is then increased by small 
increments, starting from a large negative value. After 
each increase of the field, all the spins are checked to 
find if one of them should flip (a spin flips when its ef- 
fective field changes sign). If a spin flips, its neighbors 
are checked, and so on until no spins are left that can 
flip. Then, the external field is further increased, and 
the process repeated. Since the external magnetic field 
is increased in equal increments, a large amount of time 
is spent searching the lattice for spins that can flip. The 
increments have to be big enough to avoid searching the 
lattice when there are no spins that can flip, but small 
enough so that two or more spins far apart don't flip at 
the same field. This is the method used experimentally, 
but it is suited only for "that kind of" massively parallel 
computing. 

A variation on the above method, removes the search- 
ing through the lattice that is done even if there are no 
spins that can flip. It involves looking at all the spins, 
finding the next one that will flip and then increasing the 
external field so that it does. The average searching time 
for a flip is decreased, but is still very large. Far from the 
critical point, where spins will tend to flip independently 
of each other, the time for searching scales like iV 2 where 
N is the number of spins in the system. 

The search time can be further decreased if the random 
fields are initially ordered in a list. The first spin that 
will flip is the one on "top" of the list. The external 
field is increased until the effective field of the top spin 
become zero, and the spin flips. We then check its nearest 
neighbors, and so on, while keeping the external field 
constant. When no spins are left to flip, the external field 
needs to be increased again. The change in the external 
field Aif, necessary to flip the next spin, is found by 
looking for the spin whose random field hi satisfies: 



hi > -{H Q i d + A.H") - (2n T - z)J 



(4) 



where H Q id is the field at which the previous spins have 
flipped, z is the coordination number, and nj is the num- 
ber of nearest neighbors pointing up (sj = +1) for spin 
Sj . In general, there will be a minimum of z + 1 spins to 
check from the list, since ri| can have the integer value 
between zero and z. The spin for which equation (|j) is 
satisfied for the smallest AiJ, and for which the number 
of up neighbors is n-f, will flip. In general, more than 
z + 1 spins will need to be checked because a spin can 
satisfy equation (^) for some value of rif but might not 
have that number of up neighbors, or the spin might have 
already flipped. This algorithm decreases the searching 
time since not all the spins need to be checked to find 
the next spin that will flip. Our early simulation work 
|]l3| , ^5| used this method. In practice, about half of the 
time was spent for the N log2N initial sorting of the list 
of random field numbers, where N is the total number 
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of spins in the system. The big drawback of this method 
(as for the ones mentioned above) is the huge amount 
of storage space needed to store the random fields, the 
positions of each spin, and the values of the spins. This 
becomes particularly important when larger size systems 
are simulated. 

The results in this paper use a more sophisticated algo- 
rithm which removes the need for a large storage space. 
It revolves around the idea that the change AH in the 
external field, between two avalanches, follows a proba- 
bility distribution since the random fields hi are given 
by a Gaussian distribution. The increments AH in the 
external field should be chosen according to that distri- 
bution. The probability distribution itself is not known 
explicitly, but its integral from to some finite AH is. 
It is the probability, P^° ne (AH), that no spin will flip in 
the whole system during a field change less than AH. It 
is given by: 



PT e (^H) = n„ T p:^{ah) 



(5) 



where the product is over n-f = 0, 1, z, and P™° ne (AH) 
is the probability for a down spin with nj up nearest 
neighbors not to flip when the external field changes by 
less than AH: 



(i-£ 



fl?ocal("T) 



P™ ne (AH) 



p(f) df 



P™ fUp [H local 



(6) 



The function p(f) is the random field distribution func- 
tion, and Hi oca i{n^) and HJ^^n^) are defined respec- 
tively as: 



Hi OC ai{n^) = -H - (2n T - z)J 



and 



H lo e Zl{ n ]) 



(H + AH) - (2n T -z)J. 



(8) 



P£° hp (H T local) gives the probability that a spin with rij 
up nearest neighbors has not flipped before the field has 
reached the external magnetic field value H: 



, v. 1 pH local (n T ) 

>™ fHp (H local (ni)) = 2 + J a d/ ' ( 9 ) 



P 



and N nr is the number of down spins that have n.| up 
neighbors. 

A field increment AH that has the required probabil- 
ity distribution is found by choosing a uniform random 
number between zero and one and solving for AH from 
equation (g), by setting the probability P™ ne {AH) equal 
to the value of the random number. Once the increment 
AH is known, we can find the next spin that will flip. We 
first calculate Jl^] the probability pf hp (n^) for a down 
spin with n-\ up neighbors to flip at the new field H+AH: 



P fHp (n^) 



R ni 
Rtot 



where 



P 



f P (C(«t)) 



(10) 



(11) 



is the rate at which down spins with rij up neighbors 
would flip, and Rtot is the sum of the rates R n] . for all 
nt. The spin that flips will have k up neighbors, which 
is found by satisfying the following inequality: 



£* T =0 P fUp {n ] )>O^-l P**(n T ) 



(12) 



where the cutoff C is a random number between and 
1. Once k is known, a spin is then randomly picked from 
the list of down spins with k up neighbors. 

After the first spin has flipped, its neighbors are 
checked. The probability for one of the neighbors, with 
(rij + 1) up nearest neighbors, to flip at H + AH, given 
that it has not yet flipped, is: 



P nex t(n h H + AH) = 1 



- + f 

2 ' JO 



ffiTcaiOn+i) 
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P(.f) df 



(13) 



When all the neighbors have been checked, the size of the 
avalanche is stored, as well as all the other measurements. 
The external magnetic field H is then incremented again 
by finding the next AH, starting back with equation (|^). 

The important characteristic of this method is that the 
random fields are not assigned to the spins at the begin- 
ning of the simulation, which for large system sizes de- 
creases memory requirements tremendously (asymptoti- 
cally, we use one bit per spin). This method has allowed 
us to simulate system sizes of up to 30000 2 , 1000 3 , 80 4 , 
and 50 5 spins. The majority of the data analysis was 
performed on systems of sizes 7000 2 , 320 3 , 80 4 , and 30 5 . 
The SP1 and SP2 supercomputers at the Cornell The- 
ory Center, and IBM RS6000 model 560 and J30 work- 
stations were used for the simulation. Using this new 
algorithm, close to the critical disorder, one run (for a 
particular random field configuration) for a 320 3 system 
took more than 1 CPU hour on a SP1 node at the Cornell 
Theory Center, while it took close to 37 CPU hours for a 
800 3 system on an IBM RS6000 model 560 workstation. 
Far above the critical disorder R c , the simulation time 
increases substantially: 40% above the critical disorder, 
for the 320 3 system, the simulation time was six times 
longer than for the simulation at 10% above R c . 



IV. THE SIMULATION RESULTS 

The following measurements were obtained from the 
simulation as a function of disorder R: 
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• the magnetization M(H, R) as a function of the 
external field H. 

• the avalanche size distribution integrated over the 
field H: D int (S,R). 

• the avalanche correlation function integrated over the 
field H: G mt (x,R). 

• the number of spanning avalanches N(L, R) as a 
function of the system length L, integrated over the 
field H. 

• the discontinuity in the magnetization AM(L, R) as 
a function of the system length L. 

• the second (S 2 ) rnt (L, R), third (S 3 ) in t(L, R), and 
fourth (S )i n t(L, R) moments of the avalanche size 
distribution as a function of the system length L, 
integrated over the field H. 

In addition, we have measured: 

• the avalanche size distribution D(S, H, R) as a 
function of the field H and disorder R. 

• the distribution of avalanche times Df nt ^ (S, t) as a 
function of the avalanche size S, at R — R c , integrated 
over the field H. 

The data obtained from the simulation was used to find 
and describe the critical transition. It was analyzed using 
scaling collapses. The mean field calculation for 
our model shows that near the critical point, the magne- 
tization curve has the scaling form 

M(H, R) - M C (H C , R c ) ~ \rf M ± (h/\rf s ) (14) 

where M c is the critical magnetization (the magneti- 
zation at H c , for R = R c ). r = (R c — R)/R and 
h = (H — H c ) are the reduced disorder and reduced field 
respectively, and M± is a universal scaling function (± 
refers to the sign of r). Both r and h are small. The crit- 
ical exponent (3 gives the scaling for the magnetization 
at the critical field H c (h = 0). Its mean field value is 
1/2, and the mean field value of f35 is 3/2. (Appendix A 
gives a short review on why scaling and scaling functions 
occur near a critical point, and why they have the form 
they do). 

The significance of scaling for experimental and nu- 
merical data is as follows jl8| . If the magnetization data, 
for example, is plotted against the field H, there would 
be one data curve for each disorder R (fig. ||a). While 
if we plot \r\~P M(H, R) against h/\r\^ s , all the curves 
close to R c and H c will collapse (fig. ^b) onto either 
one of two curves: one for r < (A4_), and one for 
r > (M + ). The functions M± depend only on the 
combination h/\r\P s and not on the field H and disorder 
R separately, and are therefore universal. Usually, the 
exponents are unknown and scaling or data collapses are 
used to obtain them: the exponents are varied until all 
the curves lie on top of each other. This method is useful 
for analyzing numerical as well as experimental data, and 
is often preferred to "data fitting" , as we will show. 

Numerical simulations and experiments are done on fi- 
nite size systems. Often the properties of the system will 
depend on the linear size L. Functions that depend on 



the system's length are analyzed using finite size col- 
lapses ]l8|,|l9| • An example is the number N of spanning 
avalanches: N(L,R) ~ L e N{L X / U \r\) (to be explained 
later). If N is plotted against R, there would be one data 
curve for each length L. The exponents 9 and v are ob- 
tained by plotting L~ e N(L, R) against L x l v \r\ onto one 
curve (the collapse), and extracting the exponents. 

The scaling forms we use for the collapses do not in- 
clude corrections that are present when the data is not 
taken in the limit R — > R c and L — > oo (see appendix A 
for corrections that exist in those limits). On the other 
hand, finite size effects close to R c become important. It 
is thus necessary to extrapolate to R — ► R c and L — ► oo to 
obtain the correct exponents. We have done a mean field 
simulation to test our extrapolation method. The mean 
field exponents can be calculated analytically Jb|[l4| , but 
it is useful to check that the numerical results from the 
mean field simulation, for disorders away from R c and for 
finite sizes, extrapolate to the analytical values at R = R c 
and 1/L = 0. We will see that this indeed occurs, and 
we will use the same extrapolation method in 3, 4, and 
5 dimensions. 

The mean field simulation was done with the same 
code, but with some changes. In mean field, the in- 
teractions between spins are infinite in range (each spin 
interacts when every spin in the system with the same in- 
teraction). This means that distances and positions are 
not relevant, and therefore we don't need to keep track of 
the spins and their neighbors; we just need to know the 
total number of flipped spins, and the value of the exter- 
nal field H. The following section will show the results 
of the mean field simulation and explain the extrapola- 
tion method. Then, we will turn to results in 3, 4, and 
5 dimensions. And finally, we will cover the more subtle 
scaling behavior in two dimensions. 

A. Mean Field Simulation 

The mean field simulation shows how well the results 
for the critical exponents, obtained close to R c and for 
finite size systems (finite number of spins), extrapolate 
to the calculated values for a system in the thermody- 
namic limit, at the critical disorder. Thus, we will omit 
in this section some details that are only relevant for un- 
derstanding the non-mean field simulation results. We 
start with the curves for the magnetization as a function 
of the field for different values of the disorder, which we 
find are not useful for extracting critical exponents. We 
then go on to measurements of spin avalanche sizes and 
their moments. Avalanches that span the system from 
one "side" to another will also be mentioned although 
since in mean field there are no "sides" , we will define 
what we mean by a mean field spanning avalanche. Since 
distances are irrelevant in mean field, we do not have any 
correlation measurements, but we can still apply what we 
learn from other collapses in mean field to the correlation 
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measurement data in 2, 3, 4, and 5 dimensions. 

Figure shows the magnetization curves, and figure ||a 
shows a scaling collapse for a 10 6 mean field spin system 
and r < (i? > R c ). As mentioned earlier, near the 
critical point (R c = \Jlj~K for J = 1, in mean field), the 
magnetization scales like [f3|JTT 



M{H,R) - M C (H C ,R C ) ~ \rf M ± (h/\rf s ) (15) 

where ± refers to the sign of the reduced disorder r = 
(R c — R)/R, and h = (H — H c ). The mean field critical 
exponents are (3 = 1/2 and (35 — 3/2. Notice in figure 
^a that the scaling region around M c — and H c = 
is very small; figure Bp shows that a substantially dif- 
ferent set of critical exponents leads to a similar looking 
collapse. In general, the critical field H c and the criti- 
cal magnetization M c are not zero as in mean field, and 
M c is not well determined numerically. In dimensions 
that we simulate (2 through 5), the critical region is not 
only small but it is also poorly defined, which does not 
sufficiently constrain the values of the exponents. This 
makes the magnetization function M(H, R) a poor choice 
for extracting critical exponents. 
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0.0 2.0 
Magnetic Field (H/J) 
FIG. 1. Mean field magnetization curves for 10 6 spins. 
The critical disorder is R c = 0.79788456. The curves are 
averages of 6 to 10 different initial realizations of the random 
field distribution. 

The critical magnetization M c can be removed from 
the scaling form if we look at the first derivative of the 
magnetization with respect to the field instead. dM/dH 
scales like: 



\rf-? s M ± (h/\rf s ) 



(16) 



where M± denotes the derivative of the scaling function 
Ai± with respect to its argument h/\r\^ s . The dM/dH 
mean field curves and collapses are shown in figure || 
and figures [|(a-b). Notice that the incorrect exponents 
(3 = 0.4 and (35 = 1.65 give a better collapse (fig. |]b). 
Figure || shows a close up of figure [|a, alongside with 
three (thin) curves for disorders: 0.80,0.81, and 0.82. 
These are not measured in the simulation (the finite num- 
ber of mean field spins we use give rise to finite size effects 



near R c as we will see soon) ; instead they are numerically 
calculated from the mean field implicit equation for the 
magnetization 13 14]]: 



M{H) = 1 



r -J'M(H)-H 

2 I p(f) df (17) 



where J* denotes the coupling of one spin to all the other 
spins in the system, and p(/) is the random field distri- 
bution function. 
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FIG. 2. (a) Scaling collapse for the mean field mag- 
netization curves at disorders R — 0.912, 0.974, 1.069, 
1.165, 1.197, and 1.460. (These values of disorder were chosen 
relative to R c = 0.79788456, to match some of the values we 
measured in 3 dimensions (see figure |l^). The value of the 
critical disorder R c in 3 dimensions has since been modified, 
and there is no correspondence anymore.) The exponents are 
(3 = 1/2 and (35 = 3/2. m is defined as M — M c , and in mean 
field both M c and H c are zero. The inset shows a closeup of 
the critical region, (b) Scaling collapse of the same curves as 
in (a) but with the (wrong) exponents (3 = 0.4 and (35 — 1.65. 
The two collapses are very similar. The inset is a closeup. 

The scaling collapse converges to the expected scaling 
function (dashed thick line) as we get closer to the critical 
disorder. The expected scaling form is also obtained from 
an analytic expression derived in mean field (l^Jl^l • It is 
given by the smallest real root g{y) of the cubic equation: 
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Tr 3 / 2 i? r 



y = o. 



(18) 
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0.0 

Magnetic field (H/J) 
FIG. 3. Mean field dM/dH curves for 10 6 spins and 
disorders R = 0.912, 0.974, and 1.069 (from largest to smallest 
peak). The original data is the same as in figure |l| The 
critical disorder is R c = 0.79788456. 

We again find that the critical exponents and R c , ob- 
tained from the dM/dH curves, are ill-determined. In 
finite dimensions, that is even more true since we have 
another parameter to fit: H c . For dimensions 3, 4, and 
5, we extract [3, @6, H c , and R c by other means and sim- 
ply show the resulting collapse of the M{H) and dM/dH 
curves as a check. 

As mentioned earlier, the spins flip in avalanches of 
varying sizes. The distribution of all the avalanches that 
occur at a disorders R while the external field H is raised 
adiabatically from — oo to +oo is plotted in figure §. The 
curves in this plot are normalized by the number of spins 
in the system, and therefore represent the probability per 
spin for an avalanche of size S to occur in the hysteresis 
loop, at disorder R. The curves can be normalized to one 
if they are divided by the total number of avalanches in 
the loop, and multiplied by the number of spins in the 
system. 

Often in experiments, the binned avalanche size dis- 
tribution, which contains only avalanches that occur in 
a small range of fields around a particular value of the 
field H, is measured instead. The scaling form for this 
distribution §(| is H& 
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0.0 

h/|rf 

FIG. 4. (a) Scaling collapse of mean field dM/dH 
curves from figure ^. The exponents are f3 — 1/2 and 
f35 = 3/2 (mean field values). The curves are smoothed over 
5 data points (using a running average) to show the collapse 
better. The curves do not collapse well for large and small 
h/r^ s , unless we get very close to the critical disorder (see fig- 
ure ^). (b) Scaling collapse of data in (a) but with exponents 
(3 = 0.4 and (35 = 1.65. The collapse is better, although the 
exponents are wrong. 



With the change of variable u 
becomes: 



h/\r\ ps , equation <M 



D int (S, R) ~ S 



|/35 



V ± (S"\r\,u) du (21) 



D(S,R,H) ~ S~ T V±(S a \r\ t h/\: 



(19) 



The integral in equation ( |2l| ) is a function of S a \r\ onlyso 
we can write it as: 



where S is the size of the avalanche and is large, and r 
and h are small. In mean field, a — 1/2 and r = 3/2. The 
scaling form for the integrated avalanche size distribution 
is obtained by integrating the above form over all fields: 



D int {S,R)~ / S- T V ± (S"\r\,h/\r 



dh 



(20) 



x> ( l nt) (S a \r\) 



(22) 



to obtain the scaling form for the integrated avalanche 
size distribution: 



D mt (S,R) - S~ (T+ffW X>£ nt) (S a |r|) 



(23) 



To obtain equation (22), we have assumed that the in- 
tegral in ( |2l"| ) converges. This is usually safe to do since 
the distribution curves near the critical point drop off 
exponentially for large arguments. The same kind of ar- 
gument can be made for the integrals of other measure- 
ments as well. 



G 



-3.0 -1.0 1.0 3.0 

h/|r| ps 

FIG. 5. Close-up of the mean field dM/dH curves 
collapse in figure Wa. Also plotted are three curves (thin 
lines) calculated using the mean field analytic solution to 
M(H) (see text). These are for R = 0.80, 0.81, and 0.82. 
We see that the scaling collapse, at the mean field exponents, 
of the dM/dH curves converges to the expected mean field 
scaling function (thick dashed line), as R — > R c . 




Avalanche size (S) 
FIG. 6. Mean field integrated avalanche size distri- 
bution curves for 10 6 spins and disorders R = 0.912, 0.974, 
1.069, 1.197, and 1.460 (from right to left). The straight 
line is the slope of the power law behavior in mean field: 
t + a/35 = 9/4. 

Figures and fj]b show two collapses with different 
critical exponents of the curves from figure using the 
scaling form in equation (^3|). Notice that the collapse 
with the incorrect exponents r + a (35 = 2.4 and a = 0.44 
is better than the collapse with the mean field exponents 
t + a (35 = 9/4 and a = 1/2. Although the distribution 
curves in figures [7|a and ^b have disorders that are far 
from the critical disorder R c — 0.79788456, the curves 
collapse but with the wrong exponents. 
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FIG. 7. (a) Scaling collapse of three integrated 
avalanche size distribution curves in mean held, for 

disorders: 1.069, 1.197, and 1.460. The curves are smoothed 
over 5 data points before they are collapsed. The collapse 
is done using the mean field values of the exponents a and 
r + a(3S (1/2 and 9/4 respectively), and r = (R C ~R)/R. (b) 
Same curves and scaling form as in (a), but with the expo- 
nents a — 0.44 and r + a/35 = 2.4. The collapse is better 
for the incorrect exponents! We use this "best" collapse to 
extract exponents for figures ^a and^|b, and then extrapolate 
to R = R c to obtain the correct mean field exponents. 

It is surprising that these curves collapse at all since 
the scaling form is correct only for R close to R c . Cor- 
rections to scaling become important away from the crit- 
ical point, but it seems that the scaling form has enough 
"freedom" that collapses are possible even far from R c . 
In the limit of R — > R c , we expect that the exponents 
obtained from such collapses will converge to the mean 
field value, and that the extrapolation will remove the 
question of scaling corrections. To test this, we have col- 
lapsed three curves at a time, and plotted the values of 
the exponents extracted from such collapses against the 
average of the reduced disorder \r\ for the three curves, 
which we call \r\ avg (figures ||a and |^b). In these fig- 
ures, notice two things. First, the linear extrapolation 
to \r\avg = agrees quite well with the mean field ex- 
ponent values, and second, the points obtained by doing 
collapses using r = (R c — R)/R either converge faster to 
the mean field exponents or do as well as the points ob- 
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tained from collapses done with r = (R c — R)/R c . This 
is true for all the extrapolations that we have done in 
mean field. Other models (see for example exhibit 
this behavior, and experimentalists seem to have known 
about this for a while 11231 - 
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FIG. 8. (a) r + a(3S from collapses of mean field inte- 
grated avalanche size distribution curves for 10 6 spins. 



The two points closest to 



= are for a system of 10 



spins. |r| „ 9 is the average reduced disorder \r\ for the three 
curves collapsed together (see text), (b) a from collapses of 
integrated avalanche size distribution curves for the system 
in (a). Again, the two closest points to \r\ avg = are for 
a system of 10 7 spins. The mean field values are calculated 
analytically. 

In dimensions 2 to 5, we obtain the exponents r + 
a (36 and a in the limit R — > R c , using the above linear 
extrapolation method. For other collapses, if the two 
extrapolation results differ substantially, we "bias" our 
result towards the r = (R c — R) /R extrapolated value of 
the exponent. 

Notice in figures 0a and 0b that the scaling function 
f> { l nt) has a "bump" (the — sign indicates that the col- 
lapse is for curves with R > R c ). Although we will come 
back to this point when we talk about the results in 5 
and lower dimensions, it is interesting to know what the 
shape of the scaling function p^™*- 1 is. In appendix B, 
we calculate the mean field scaling function for r < 



(equations B14 and B15| ) 



V { l nt \-rS a ) = ' ' 



(-(-r)S° u-4) 



TT-v/2 



du 



(-rS a + u) ^ 



(24) 
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FIG. 9. Scaling collapse of the mean field integrated 
avalanche size distribution curves (dashed lines), for 
S = 10 6 spins and R = 0.912,0.974, and S = 10 7 spins 
and R = 0.854,0.878,0.912. The critical exponents are: 
r + a (38 = 9/4 and a = 1/2. The thick black line is the 
best fit to the data using a function that is the product of a 
polynomial and an exponential (eqn. (pB[)). The thick grey 
line is the "real" mean field scaling function (see text). 

A closed analytic form can not be obtained, but we can 
find the behavior of this function for small and large argu- 
ments —rS 17 . For small arguments X = —rS' 7 , the scal- 



ing function is a polynomial in X ( B17 ), while for large 
arguments, the scaling function is given by the product 
of an exponential decay in X 2 and the square root of X 
(B19). We can then try to fit our data (the scaling col- 
lapse) with a function that will incorporate a polynomial 
and an exponential decay (as an approximation to the 
real function). We obtain: 



(0.204 + 0.482AT - 0.391AT 2 
0.204X 3 - 0.048X 4 ) 



(25) 



This form has the expected exponential behavior at large 
X, but the wrong pre-factor. On the other hand, for 
small X, the above function is analytic. A better ap- 
proach might be to use a parametric representation p3j , 
which we have not yet tried. 

Equation (^5|) can be compared with the curve ob- 
tained by numerically integrating the scaling function 
V { i nt) in equation @. Fi gure 9 shows the fit in black 
(equation (p^)) to the collapsed data, for curves (dashed 
lines) of different disorder, and system size S — 10 6 and 
S = 10 7 spins. The grey curve is the "real" scaling func- 
tion obtained from the numerical integration of equation 
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(p4|). Notice that the scaling collapse (done with the 
mean field values of the exponents T + af3S and a) of even 
a system of 10 7 spins and within 7% of R c (ie. R — 0.854) 
(this is the curve with the smallest peak in the graph) is 
not close to the "real" scaling function (the thick grey 
curve). The error is within 5% for this curve (within 
10% for the fit). However, as R — ► R c , the avalanche size 
distribution curves seem to be approaching the "real" 
scaling function (grey curve). It is important to keep 
in mind when analyzing experimental or numerical data 
as we will in 5 and lower dimensions, that the scaling 
collapse most likely does not give the limiting curve one 
would obtain for 1/L — ► and R — > R c , even for what 
seems like a large size, and close to the critical disorder. 
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FIG. 10. Number of mean field spanning avalanches 

N m f as a function of the disorder R. Curves at sizes 125 
and 343 are not plotted. All the error bars are not shown for 
clarity. The ones that are shown are representative for the 
peaks. The error bars are smaller for larger disorders. About 
26 points are used for each curve; each point being an average 
between 250 (for size 512000) and 2500 (for size 1000) random 
field configurations. The inset shows the collapse of the three 
largest size curves using the mean field (calculated) exponents 
6 = 3/8 and 1/v = 1/4. 

The avalanches in the avalanche size distribution are 
finite, by which we mean that they don't span the sys- 
tem. We have mentioned earlier that due to the fi- 
nite size of a system, close to the critical disorder R c , 
the largest avalanche or avalanches will span the system 
from one side to another. We will talk about spanning 
avalanches in more details later, but for now we just need 
to know that the number N of spanning avalanches scales 
as N(L, R) ~ L e N±(L 1 / V \r\) where N± is a scaling func- 
tion (± indicates the sign of r), L is the linear size of the 
system, 9 is the exponent that arises from the existence 
of more than one spanning avalanche, and v is the corre- 
lation length (£) exponent: £ ~ \r \~ p . 
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FIG. 11. (a) and (b) 1/u and 6 respectively, extracted 
from the mean field spanning avalanches collapses, as 

a function of the geometric average of 1/ S m f for three curves 
collapsed together (see text). The extrapolation (non-linear 
for 8) to 1/S m f ™ > agrees with the calculated values for the 
two exponents. 

As was mentioned earlier, in mean field there is no 
meaning to distance or lattice, and thus there are no 
"sides" . Purely for the purpose of testing our extrapola- 
tion method for finite size scaling collapses in the mean 
field simulation, we have defined a mean field "spanning 
avalanche" to be one with more than \J S m f spins flip- 
ping at a field H , where S m t is the total number of spins 
in the system. (Note that the mean field exponents are 
valid for dimensions 6 and above, but that in those di- 
mensions distances do have a meaning.) Using the above 
definition of a mean field spanning avalanche, it can be 
shown (see appendix C) that the scaling form for their 
number is: 



N mf (S mf ,R)~S« mf N? f (S 1 J;\r 



(26) 



and that the values of the exponents 6 and 1/u are 3/8 
and 1/4 respectively. N m f is the number of mean field 
spanning avalanches, while A/J 1 is a universal scaling 
function. The exponents 9 and 1/v are defined by the 
arbitrary definition for a spanning avalanche. Because 
of how they are defined, their values are different from 
the mean field values of 1/v = 2 and 6 = 1, obtained 
from the renormalization group [|lj|JTq] and the exponent 
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scaling relation 1/er = (d — 9)v — (3 [^4 16 1. 

Figure [To] shows the number of mean field spanning 
avalanches as a function of disorder, for several sizes, 
as well as the scaling collapse of the data. Note that 
the number of spanning avalanches close to the criti- 
cal disorder R c = increases with the size S m f 
of the system, and that the peaks are getting narrower. 
The scaling collapse in the inset, shows only the three 
largest curves. For smaller sizes, the peaks do not col- 
lapse well with the larger size systems presumably due 
to finite size effects. The extrapolation plots for 9 and 
1/9 are shown in figures [Tl| a and pd|b. On the horizontal 
axis of these two plots is the geometric mean of 1/ S m f 
for the three curves that are collapsed together, analo- 
gous to the extrapolation method used for the integrated 
avalanche size distribution. Note that the extrapolation 
to l/S m f — > for 8 does not seem to be linear, and 
that the value of \jv from the linear extrapolation of the 
r = (R c — R)/R data agrees better with the mean field 
value than the value obtained from the linear extrapola- 
tion of the r — {R c — R) /R c data. 

Note that we measure the avalanche size distribution 
only for disorders at which there are no "mean field span- 
ning avalanches" (for a 10 6 system, that is for R > 0.912), 
since that's what we do in dimensions 2 through 5 (finite 
dimensions) to avoid large finite size effects. For the sec- 
ond moments of the avalanche size distribution measure- 
ments (see below) , the spanning avalanches were removed 
(same as in finite dimensions). 

We have also measured the change in the magnetiza- 
tion AM due to all the spanning avalanches, as a function 
of the disorder R (figure p"2|a). This gives us an indepen- 
dent measurement of the exponent (3. In the thermo- 
dynamic limit above the critical disorder, there are no 
spanning avalanches so the change in the magnetization 
AM will be zero, while for small disorders the change 
in the magnetization will converge to one. Close and 
below the critical disorder R c , at the critical field, the 
scaling form for the change in the magnetization due to 
the spanning avalanches will be (from equation (|15|) 



AM(H = H C ,R) 



\0 



(27) 



For finite size systems, as shown in the figure, the change 
in the magnetization is not zero above the critical disor- 
der: the data has to be analyzed using finite size scaling. 
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FIG. 12. (a) Change in the magnetization due to 

spanning avalanches as a function of disorder R. The data 
is for several mean field system sizes. The critical disor- 
der is R c = 0.79788456. The statistical errors are not 
larger than 0.005 (in units of AM), (b) Mean field scal- 
ing collapse of the change in the magnetization curves for 
sizes S m f = 1000,8000,64000,512000. The exponents are 
l/u = 0.25 and (3/v = 0.125 and r = (R c - R)/R. The part 
of the curve that is collapsed is for R> R c . 

The dependence on the system size S m f can be brought 
in through a scaling function (see references jl8],[l9| ) that 
we call AA4±: 



AM(S mf ,R) 



AM±(S, 



1/5, 
mf I 



(28) 



where v is defined above, and ± refers to the sign of r. 
We are free to define the scaling function AM.± as: 



AM ± (S^\r\) ee (Ol) " AM ± (S^\r\), (29) 

where AA4± is now a different scaling function. The 
scaling form for the change of the magnetization AM 
then becomes: 



AM(S mf ,R) ~ S-V* AM±(S^\r\) 



(30) 



Figure |12| b shows a collapse of the data using this scaling 
form. The collapse is done for disorders close to and 
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above the critical disorder, that is, for r < 0. The scaling 
function in figure |i"^b, in the range of the collapse, is 
therefore AA4_. 
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FIG. 13. (a) and (b) Mean field exponents \jv and f3/v 
respectively, from collapses of the magnetization change 
due to spanning avalanches (see text). The extrapolation to 
(1/S m f)gm = agrees with the calculated values. 

Values for the exponents 1/v and (3/v extracted from 
such collapses at several geometric average reciprocal 
sizes are shown in figures fl3| a and |i~3|b. (These plots are 
done the same way as for the spanning avalanches expo- 
nents.) The linear extrapolation to 1/S m f = is in very 
good agreement with the calculated values. Note that 
the extrapolation for l/i> of the r = (R c — R)/R data 
gives again a better agreement with the calculated value 
than the extrapolation using the r — (R c — R) /R c data. 
The exponent (3 in 3, 4, and 5 dimensions is calculated 
from (3/v, which is extracted from the above kind of col- 
lapse. The obtained value is used to check the collapse 
of the M(H) and dM /dH data curves. 
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FIG. 14. Mean field second moments of the avalanche 
size distribution integrated over the field H, for several differ- 
ent sizes. More than 20 points are used for each curve; each 
point being an average of a few to several hundred random 
field configurations. The error bars for the S m f ~ 1, 000 curve 
are too small to be shown. Curves at S m f = 125 and 343 are 
not shown. The inset shows the collapse of these four curves 
at l = — (t + crPS — 3)/ai> = 3/8 and 1/v = 1/4, which are 
the mean field calculated values. 

Another quantity that is related to the avalanches is 
the moment of the size distribution. We have measured 
the second, third, and fourth moment, and we will show 
how the second moment scales and collapses in mean 
field. The second moment is defined as: 



(S 2 



S 2 D{S, R, H, S, 



mf 



dS 



(31) 



where D(S, R, H, S m f) is the avalanche size distribution 
mentioned above, but with the system size S m f included 
as a variable since we are looking for the finite size scal- 
ing form, as is clear from the data in figure [l4|. Recall 
that only non-spanning avalanches are included in the 
distribution function D(S, R, H, S m f). Equation |^ can 
be written in terms of the scaling form for large sizes S 
of the avalanche size distribution D: 



(S 2 



S 2 S~ T V ± (S a \r\,h/\; 



S]S\r\ 



dS (32) 



As we have seen before, the dependence on the system 
size in the scaling function f>± is given by S^J\r\ where 
v is defined above through the definition of a mean field 
spanning avalanche. If we define 



p ± (s>i,vkr<f>i) = 

(^|r|)^ V ± (S^\r\,h/\rf s ,S^\r\) 



(33) 



where T>± is a different scaling function, and let u = 
Sir] 1 /*, we obtain: 

(S 2 ) ~ \r\^-V/° Jv ± (u\h/\rf s ,S^\r\)du (34) 

The integral in equation ( pi[ ) is a function of h/\r\P 5 and 
Smf l r l on ly> so we can wr he: 
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(S 2 



,(r-3)/ CT c(2) 



S£>(h/\rf s ,S%\r\) 



— (2) 

(35) where S± is a universal scaling function (± indicates the 



(2) 

which is the second moment scaling form, and S± is a 
universal scaling function. 
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FIG. 15. (a) and (b) Values for and (r + a/35-3)/(ri> 
respectively, extracted from the collapses of the second mo- 
ments of the avalanche size distribution. The exponents are 
plotted as a function of the geometric average of 1/ S m f for 
three curves collapsed at a time (see text). The extrapola- 
tion to large sizes agrees with the calculated values for these 
exponents. 

In the simulation, we have measured the second mo- 
ment of the distribution integrated over the field H, 
whose scaling form can be obtained by integrating the 
result of equation (|35"|): 

(S 2 ) mt ~ M (r - 3)/CT | s¥\h/\rf s ,Sl(;\r\) dh (36) 

As was done previously, we define u = h/\r\@ s , and call 
the remaining integral: 

Sg\h/\rf s ,S$\r\)dh = 
(S^\r\)-^ S - 3 ^ S^\S^\T\) (37) 

to obtain the second moment of the avalanche size dis- 
tribution integrated over the magnetic field H: 



-(T+cr/3(5-3)/cr£i ^(2) (Q l/v 
int ~ u mf 



(38) 



sign of r). The mean field value for — (r + cr(3S 
is 3/8. 



3)/cri> 



Figure 14 shows the integrated second moments of 



non-spanning mean field avalanches for several system 
sizes, and a collapse using the scaling form in equation 
(|38|). Figures p"5| a and [T^ b show the values for l/i> and 
— (t + a/35 — i)/av respectively, for several geometric av- 
erage reciprocal sizes, and show how well they linearly 
extrapolate to l/S m f — ► 0. These plots are done the 
same way as for the mean field spanning avalanches. No- 
tice that for l/z>, the linear extrapolation of the data 
using r = (R c — R) / R gives a much better agreement 
with the calculated value than the linear extrapolation 
of the data obtained using r = (R c — R)/R c . 

To summarize this section, we have shown that the 
values of the critical exponents extracted from our mean 
field simulation by scaling collapses, extrapolate to the 
expected (calculated) values for R — > R c and 1/S m f — > 0. 
Thus corrections to scaling due to finite sizes as well as fi- 
nite size effects near the critical point seem to be avoided 
by extrapolation. The same extrapolation method is 
therefore used for extracting exponents in 3, 4, and 5 
dimensions, which we will see next. The results in 2 di- 
mensions will be shown last. 



B. Simulation Results in 3, 4, and 5 Dimensions 

1. Magnetization Curves 

The magnetization as a function of the external field 
H is measured for different values of the disorder R. Ini- 
tially all the spins are pointing down (sj = — 1 for all 
i). The field is then slowly raised from a large negative 
value, until a spin flips. When the first spin has flipped, 
the external field is held constant while all the spins in 
the avalanche are flipping. The change in the magneti- 
zation due to this avalanche is just twice the size of the 
avalanche. 

Figure fl6| a shows the magnetization curves obtained 
from our simulation in 3 dimensions for several values of 
the disorder R. Similar plots can be obtained in 4 and 
5 dimensions p^j . As the disorder R is decreased, a dis- 
continuity ("jump") in the magnetization curve appears. 
The critical disorder R c is the value of the disorder at 
which this discontinuity appears for the first time as the 
amount of disorder is decreased, for a system in the ther- 
modynamic limit. For finite size systems, like the ones 
we use in our simulation, the "jump" will occur earlier. 
The effective critical disorder for a system of size L is 
larger than the critical disorder R c (1/L = 0). The criti- 
cal disorder R c is found from finite size scaling collapses 
of the spanning avalanches and second moments of the 
avalanche size distribution which will be covered later. 
The values are listed in Table ffl. 
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We have seen in mean field that the magnetization 
curves near the transition scale as 



1.0 



m(H,R) ~ \rf M ± (h/\rf d ) 



(39) 



where m = M(H,R) - M C (H C ,R C ), h = H - H c , and 
M± is the corresponding scaling function. The critical 
magnetization M c and critical field H c are not universal 
quantities: in our mean field simulation and the hard- 
spin mean field model for our system fl4j| , both are zero; 
however they are non-zero quantities in a soft-spin model 

In general, the scaling variables in (p9l) need not be 
r and h. but can instead be some "rotated" variables r' 
and h! fl25|| which to first approximation can be written 
as: 



and: 



ah 



ti =h + br 



(40) 



(41) 



(See appendix A for these and other corrections.) The 
constants a and b are not universal and the critical ex- 
ponents do not depend on them (for the mean field data 
a = b = 0). In equation (|39"|), the scaling variables r and 
h should be replaced by the "rotated" variables r' and 
hi ' , but since the measurements in our simulation are in 
terms of r and h, we rewrite the scaling form in terms 
of those. We find that in the leading order of scaling 
behavior, the magnetization scales like: 



M(H, R) - M c 



I u 



M- 



:({h + br)/\r 



fib 



(42) 



The correction br is dominant for R — > i? c , and can not 
be ignored. The opposite is true for a h (see appendix 
A). 

From the previous equation, the parameters that need 
to be fitted are M c , H c , (3, (38, and the "tilting" constant 
b. These should be found by collapsing the magnetization 
curves onto each other. As in mean field, we find that 
collapses of magnetization curves in 3, 4, and 5 dimen- 
sions do not define well the value of the critical magneti- 
zation M c . Furthermore, we observe strong correlations 
between the parameters, which lead to weak constraints 
on their values. 
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FIG. 16. (a) Magnetization curves in 3 dimensions 

for size L = 320, and three values of disorder. The curves 
are averages of up to 48 different random field configurations. 
Note the discontinuity in the magnetization for R = 2.20. 
In finite size systems, the discontinuity in the magnetization 
curve occurs even for R > R c (R c = 2.16 in 3 dimensions), 
(b) "Tilted" scaling collapse (see text) of the magnetization 
curves in 3 dimensions for size L = 320. The disorders range 
from R = 2.35 to R = 3.20 (7? > R c ). The critical magne- 
tization is chosen as M c = 0.9 from an analysis of the mag- 
netization curves and is kept fixed during the collapse. The 
exponents are /3 = 0.036, f38 — 1.81, and the critical field 
and disorder are 1.435 and 2.16 respectively. The "tilting" 
parameter b is 0.39. 

To remove the dependence on the critical magnetiza- 
tion M c , we can look at the collapse of dM/dH which 
scales like: 



\r\P~P s M±({h + br)/\ 



IPS 



(43) 



Although M c does not appear in the above form, the 
other parameters are still not uniquely defined by the 
collapse. We find that we need to extract (3 from the 
magnetization discontinuity (AM) collapses, and (3d and 
H c from the binned avalanche size distribution collapses 
rather than from the magnetization curves themselves. 
Using the values obtained from these collapses, and the 
value of R c , the "tilting" constant b is then found from 
magnetization curve collapses (figure |l6|b) . 
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FIG. 17. (a) Derivative with respect to the field H 
of the magnetization M, for disorders R = 2.35, 2.4, 2.45, 
2.5, 2.6, 2.7, 2.85, 3.0, and 3.2 (highest to lowest peak), in 
3 dimensions. The curves are smoothed by 10 data points 
before they are collapsed, (b) Scaling collapse of the data in 
(a) with ft = 0.036, PS = 1.81, 6 = 0.39, H c = 1.435, and 
R c = 2.16. While the curves are not collapsing onto a single 
curve, neither did they for the mean field theory curves (figure 
^,). This is because the curves are still far from the critical 
disorder R c . 



Figure |17| a shows the curves for the derivative of the 
magnetization with respect to the field H, and figure [j~7| b 
shows the scaling collapse using the same exponent and 
parameter values as in figure [T6|b. The collapsed curves 
have disorders larger than the critical disorder: below 
R C) the fluctuations are larger and the collapses are less 
reliable. 

Since we found that b ^ (b = 0.39 in 3d), the scaling 
variables are indeed some r' and h! , and not the variables 
we measure: r and h. Therefore, the scaling functions 
will in general be functions of a different combination of 
scaling variables from the ones we used in mean field, 
where the scaling variables are r and h. However, we 
find in appendix A that the measurements that are inte- 
grated over the external field H remove the "tilt" param- 
eter b (other analytic corrections might still be important 
though). This is true for the integrated avalanche size 
distribution, the avalanche correlation (integrated over 
the field), the number of spanning avalanches, the mo- 



ments of the avalanche size distribution, and the time 
distribution of avalanche sizes. In the sections that treat 
these measurements, we will ignore the "rotation" of axis 
to simplify the presentation. Note that the change in 
the magnetization AM due to the spanning avalanches 
is integrated over only a small range of external fields 
(wherever there are spanning avalanches). On the other 
hand, the binned avalanche size distribution is not inte- 
grated over the field H, and we therefore examine this 
measurement more carefully. 



2. Avalanche Size Distribution 

a. Integrated Avalanche Size Distribution In our 
model the spins often flip in avalanches, which are col- 
lective flips of neighboring spins at a constant external 
field H . These avalanches come in different sizes. The 
integrated avalanche size distribution is the size distri- 
bution of all the avalanches that occur in one branch of 
the hysteresis loop (for H from — oo to oo). Figure [l8| 
p"2[ shows some of the raw data (thick lines) in 3 dimen- 
sions. Note that the curves follow a power law behavior 
over several decades. Even 50% away from criticality (at 
R = 3.2), there are still two decades of scaling, which 
implies that the critical region is large. In experiments, 
a few decades of scaling could be interpreted in terms 
of self-organized criticality (SOC). However, our model 
and simulation suggest that several decades of power law 
scaling can still be present rather far from the critical 
point (note that the size of the critical region is non- 
universal). In the figure, the cutoff in the power law 
which diverges as the critical disorder R c is approached 
(R c = 2.16 in 3 dimensions), is a signature that the sys- 
tem is away from criticality, and that a parameter can be 
tuned (here R) to bring it to the transition. This cutoff 
scales as S ~ |r|~ 1 / <T , where S is the avalanche size and 
r = (R c — R)/R is the reduced disorder. 

The power law for the curves of figure [l^ can be ob- 
tained through scaling collapses. A plot is shown in the 
inset of figure |l^. The scaling form is (see mean field 
section) 



Di n t(S, R) ~ S 



-(r+a/36) j){int) rgo 



(44) 



where l?^™*- 1 is the scaling function (the — sign indicates 
that the collapsed curves are for R > R c ). The critical 
exponents t + a/35 = 2.03 and a = 0.24 are obtained 
from collapses and linear extrapolation of the extracted 
values to R — R c (figures |H^ a and 19 b), as was done in 
mean field. (Although the "real" scaling variables are r' 
and ft', when integrating over the field H we recover the 
same form as in mean field; see appendix A.) Table || 
lists all the exponents extracted from scaling collapses, 
and extrapolated to R — ► R c and 1/L — > 0. 
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Avalanche size (S) 
FIG. 18. Avalanche size distribution integrated over 
the field H in 3 dimensions, for 320 3 spins and disor- 
ders R = 4.0, 3.2, and 2.6. The last curve is at R = 2.25, 
for a 1000 3 spin system. The 320 3 curves are averages over 
up to 16 initial random field configurations. All curves are 
smoothed by 10 data points before they are collapsed. The in- 
set shows the scaling collapse of the integrated avalanche size 
distribution curves in 3 dimensions, using r = (R c — R)/R, 
t + a (35 = 2.03, and a = 0.24, for sizes 160 3 , 320 3 , 800 3 , 
and 1000 3 , and disorders ranging from R — 2.25 to R — 3.2 
(R c = 2.16). The two top curves in the collapse, at R = 3.2, 
show noticeable corrections to scaling. The thick dark curve 
through the collapse is the fit to the data (see text). In the 
main figure, the distribution curves obtained from the fit to 
the collapsed data are plotted (thin lines) alongside the raw 
data (thick lines). The straight dashed line is the expected 
asymptotic power law behavior: S~ 2 0S , which does not agree 
with the measured slope of the raw data due to the shape of 
the scaling function (see text). 

We have mentioned earlier that the mean field scal- 
ing function f> { i nt) (X), with X = S a \r\ and r < 0, is a 
polynomial for small X and gives an exponential in X 1 ! 
(l/cr = 2 in mean field) multiplied by X 13 ((3 — 1/2 in 
mean field) for large X (see mean field section and ap- 
pendix B). As we have done in mean field, we can try to 
fit the scaling function X>^ nt ^ in dimensions 5 and below 
with a product of a polynomial and an exponential func- 
tion. This is done in 3 dimensions in the inset of figure 
[l8| (thick black line through the data). The phenomeno- 
logical fit is: 

(0.021 + 0.002X + 0.531X 2 - 0.266X 3 + 0.261X 4 ) (45) 

with l/cr = 4.20 which is obtained from scaling collapses. 
The distribution curves obtained using the above fit are 
plotted (thin lines in figure |l8|) alongside the raw data 
(thick lines) . They agree remarkably well even far above 
R c . We should recall though, from the mean field discus- 
sion (see figure [)]) , that the fitted curve to the collapsed 
data can differ from the "real" scaling function even for 
large sizes and close to the critical disorder (in mean field 
the error was about 10%). We expect a similar behavior 
in finite dimensions. 
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FIG. 19. (a) and (b) r + a/38 and a respectively, from 
collapses of the integrated avalanche size distribution 
curves for a 320 3 spin system. The data is plotted as in 
mean field. The two closest points to |r| TOg =0 are for a 800 3 
system, for a collapse using curves with disorder: 2.26, 2.28, 
2.30, 2.32, 2.34, and 2.36. The extrapolation to \r\ avg = 
gives: r + a/3S = 2.03 and a = 0.24. 

The scaling function in the inset of figure [l8] has a pe- 
culiar shape: it grows by a factor often before cutting off. 
The consequence of this shape is that in the simulations, 
it takes many decades in the size distribution for the 
slope to converge to the asymptotic power law. This can 
be seen from the comparison between a straight line fit 

and the 
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through the R — 2.25 (1000 3 !) curve in figure 
asymptotic power law S~ 2 03 obtained from scaling col- 
lapses and the extrapolation (thick dashed straight line 
in the same figure). A similar "bump" exists in other di- 
mensions and mean field as well. Figure ||(] shows the 
scaling functions in different dimensions and in mean 
field. In this graph, the scaling functions are normal- 
ized to one and the peaks are aligned (the scaling forms 
allow this) . The curves plotted in figure |o| are not raw 
data but fits to the scaling collapse in each dimensions, 
as was done in the inset of figure [l8| The mean field 
and 3 dimensions curves are given by equations ( pq ) and 
d45| ) respectively. For 5, 4, and 2 dimensions, we have 
respectively: 



V { * nt) (X) = 



-0.518X 1/o 
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(0.112 + 0.459X - 0.260X 2 + 0.201X 3 - 0.050X 4 ) (46) 



v { l nt \x) 



-0.954X 1/o 



(0.058 + 0.396X 



V [ l nt) {X) = e' imexl " x 
(0.492 - 4.472X + 14.702X 2 - 
20.936X 3 + 11.303X 4 ) 



0.248X 2 - 0.140X 3 



0.026X 4 ) (47) 



(48) 



with 1/cr = 2.35, 3.20, and 10.0. The errors in the fits are 
in the same range as for the mean field simulation data 
(see figure ^). The 2 dimensional fit plotted in grey will 
be covered further in the next section. 
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FIG. 20. Integrated avalanche size distribution scal- 
ing functions in 2, 3, 4, and 5 dimensions, and mean 
field. The curves are fits (see text) to the scaling collapses 
done with exponents from Table [n] and VII. The peaks are 
aligned to fall on (1,1). Due to the "bump" in the scaling 
function the power law exponent can not be extracted from a 
linear fit to the raw data. 

From figure ^0] we can conclude that in each dimension 
(and in mean field!), a straight line fit to the integrated 
avalanche size distribution data is going to give the wrong 
critical exponent, and that only by doing scaling collapses 
and an extrapolation the asymptotic value can be found. 
This is shown for 3 dimensions in figure ||, and was 
found to be true in other dimensions as well. We will 
next see that this is different for the binned avalanche 
size distribution. The value for the slope obtained from 
a linear fit to the data agrees very well with the value 
obtained from the scaling collapses. 
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FIG. 21. (a) Binned in H avalanche size distribution 
in 4 dimensions for a system of 80 4 spins at R = 4.09 
(Re = 4.10). The critical field is H c = 1.265. The curves are 
averages over close to 60 random field configuration. Only 
a few curves are shown, (b) Scaling collapse of the binned 
avalanche size distribution for H < H c (upper collapse) and 
H > H c (lower collapse). The critical exponents are r = 1.53 
and cr/35 = 0.54, and the critical field is H c = 1.265. The bins 
are at fields: 1.162, 1.185, 1.204, 1.220, 1.234, 1.245, 1.254, 
1.276, 1.285, 1.296, 1.310, 1.326, 1.345, and 1.368. Notice that 
the two scaling functions do not have a "bump" (see text). 

b. Binned in H Avalanche Size Distribution The 
avalanche size distribution can also be measured at a field 
H or in a small range of fields centered around H. Wc 
have measured this binned in H avalanche size distribu- 
tion for systems at the critical disorder R c (r = 0). To 
obtain the scaling form, we start from the distribution of 
avalanches at field H and disorder R (eqn. |l9|): 



D(S,R,H)~S- T V ± (S"\r\,\h\/\rf s ) 



(49) 



where as before T>± is the scaling function and ± indi- 
cates the sign of r. (For most of our data, we can ignore 
the corrections due to the "rotation" of axis as explained 
in appendix A.) The scaling function can be rewritten as 

t> ± (S a \r\, {S a \r\y s \h\/\r\^, where D± is a new scaling 
function. Letting R — > R c , the scaling for the avalanche 
size distribution at the field H, measured at the critical 
disorder R c is: 
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D{S,H) ~ S~ T V ± (\h\S a ^) 



(50) 
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FIG. 22. Values for the exponent r extracted from 
the binned in H avalanche size distribution curves in 4 
dimensions, for a 80 4 spin system at R = 4.09 (R c = 4.10). 
The critical field is H c — 1.265. The exponent r is found 
from this linear extrapolation to AH avg — 0. The exponent 
a(3& is calculated from the value of r + a/35, extracted from 
the integrated avalanche size distribution, and the value of r 
from this plot. 

Figure pl| a shows the binned in H avalanche size dis- 
tribution curves in 4 dimensions, for values of H below 
the critical field H c . (The curves and analysis are sim- 
ilar in 3 and 5 dimensions; results in 4 dimensions are 
used here for variety.) The simulation was done at the 
best estimate of the critical disorder R c (4.1 in 4 dimen- 
sions) . The binning in H is logarithmic and started from 
an approximate critical field H c obtained from the mag- 
netization curves; better estimates of H c are obtained 
from the binned distribution data curves and their col- 
lapses. Our best estimate for the critical field H c in 4 
dimensions is 1.265 ± 0.007. The scaling form for the 
logarithmically binned data is the same as in equation 
( |50| ) , if the log-binned data is normalized by the size of 
the bin. Figure pl| b shows the scaling collapse for our 
data, both below and above the critical field H c . The 
"top" collapse gives the shape of the X>_ (H < H c ) func- 
tion, while the "bottom" collapse gives the T> + (H > H c ) 
function. Above the critical field H Cl there are spanning 
avalanches in the system [^6). These are not included in 
the binned avalanche size distribution collapse shown in 
figure pl|b. 

The exponent r which gives the power law behavior of 
the binned avalanche size distribution is obtained from 
an extrapolation similar to previous ones (figure ^2|) , but 
with the field H (AH avg in figure |2^ is the algebraic 
average of H — H c for three curves collapsed together) 
as the variable instead of the disorder R. The exponent 
a (36 is found to be very sensitive to H c , while r is not. 
We have therefore used the values of r + a (3d and a from 
the integrated avalanche size distribution collapses, and 
t from the binned avalanche size distribution collapses 



to further constrain H c (by constraining a (35), and to 
calculate (38. The latter is then used to obtain collapses of 
the magnetization curves. We should mention here that 
H c in all the dimensions is difficult to find and that it is 
influenced by finite sizes. The values listed in Table [| are 
the best estimates obtained from the largest system sizes 
we have. Nevertheless, systematic errors for H c could 
be larger than the errors given in Table |. This implies 
possible systematic errors for a (35 which depends on H c , 
and for (36 which is calculated from a(36. These could 



also be larger than the errors listed in Table III 
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FIG. 23. Binned avalanche size distribution curve 
(dashed line) in 4 dimensions, for a system of 80 4 spins 
at R c = 4.09. The magnetic field is H = 1.265. The straight 
solid line is a linear fit to the data for S < 13, 000 spins. The 
slope from the fit is 1.55 (this varies by not more than 3% 
as the range over which the fit is done is changed), while the 
exponent r obtained from the collapses and the extrapolation 
in figure p2| is 1.53 ±0.08. 



From figure |21|b, we see that the two binned avalanche 
size distribution scaling function do not have a "bump" 
as does the scaling function for the integrated avalanche 
size distribution (inset in figure |l§| ). Therefore, we ex- 
pect that the exponent r which gives the slope of the 
distribution in figure pl| a can also be obtained by a lin- 
ear fit through the data curve closest to the critical field. 
Figure ^ shows the curve for the H = 1.265 bin (dashed 
curve) as well as the linear fit. The slope from the linear 
fit is 1.55 while the value of r obtained from the collapses 
and the extrapolation in figure ^ is 1.53 ± 0.08. Fitting 
the binned distribution curves with a straight line to ex- 
tract the exponent r is also possible in other dimensions 
and mean field as well. 



3. Avalanche Correlation 

The avalanche correlation function G(x, R, H) mea- 
sures the probability that a flipping spin will trigger, 
through an avalanche of spins, another spin a distance 
x away. From the renormalization group description 
0,fl5[, close to the critical point and for large distances 
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x, the correlation function is given by (corrections are 
subdominant as explained in appendix A): 



G(x,R, H) 



1 



d-2+7) 



g±(x/Z(r,h)) 



(51) 



where r and h are respectively the reduced disorder and 
field, Q± (± indicates the sign of r) is the scaling function, 
d is the dimension, £ is the correlation length, and 
is called the "anomalous dimension". The correlation 
length £ (r, h) is a macroscopic length scale in the system 
which is roughly on the order of the mean linear extent of 
the avalanches for a system away from the critical point. 



10 J 



3 10 ° 



10" 




10 10' 10' 

Distance (x) 

FIG. 24. Avalanche correlation function integrated 
over the field J? in 3 dimensions, for L — 320. The 
curves are averages of up to 19 random field configurations. 
The critical disorder R c is 2.16. 

At the critical field H c (h=0) and near R c , the correla- 
tion length scales like £ ~ \ r \~ L ' > while for small field h it 
is given by £ ~ \ r \~ u y±(h/\r\^ s ) where y± is a univer- 
sal scaling function. The avalanche correlation function 
should not be confused with the cluster or "spin-spin" 
correlation which measures the probability that two spins 
a distance x away have the same value. (The algebraic 
decay for this other, spin-spin correlation function at the 
critical point (r = and h = 0), is l/x d ~ 4+j i ftl.) 
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FIG. 25. (a) Scaling collapse of the avalanche correla- 
tion function integrated over the field H, in 3 dimen- 
sions for L — 320. The values of the disorders range from 
R — 2.35 to R — 3.0, with R c = 2.16. The exponents are: 
v = 1.39 ± 0.20 and d + j3/v = 3.07 ± 0.30. (b) Exponent v 
extracted from collapses of avalanche correlation curves (see 



(a)). The extrapolated value at 



= is 1.37 ±0.18. 



We have measured the avalanche correlation function 
integrated over the field H, for R > R c . For every 
avalanche that occurs between H = — oo and H = +oo, 
we keep a count on the number of times a distance x 
occurs in the avalanche. To decrease the computational 
time not every pair of spins is selected; instead we do 
a statistical average for S pairs where S is the size of 
the avalanche. Our simulation seems to indicate that the 
difference between this statistical average and the exact 
measurement is less than the fluctuations obtained from 
measurements of the correlation function for different re- 
alizations of the random field distribution. The data is 
saved in "distance" bins separated by 0.5 and starting 
at a distance of 1.0 (the self correlation is not included), 
and is normalized by the number of neighbors at each 
distance. The spanning avalanches are not included in 
our correlation measurement. Figure [24| shows several 
avalanche correlation curves in 3 dimensions for L = 320. 
The scaling form for the avalanche correlation function 
integrated over the field H, close to the critical point and 
for large distances x, is obtained by integrating equation 
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G mt (x,R) ~ / g±(x/£(r,h)) dh (52) 

Near the critical point £(r, |r-|^ ^3^ ± (^/ |r- 1 ■ s<5 ) . Defin- 
ing u = h/\r\@ s , equation ( p2[ ) becomes: 

G mt (x,R) ~ Irl^x"^-^ /" ^(Vkl""^^)) d« 

(53) 

The integral (X) in equation ( |53] ) is a function of x|r| w 
and can be written as: 

l=(x\r\ v )-f }S ' v Q±(x\r\ v ) (54) 

to obtain the scaling form: 

G ™( x > R )~-^j3frS±(x\r\ v ) (55) 

where we have used the scaling relation [2 — rf)y = (35 — (3 
(see 0,^6) for the derivation). 
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FIG. 26. Anisotropies in the avalanche correlation 
function. The curves are for a system of 320 3 spins at 
R — 2.35. Four curves are shown on the graph: one is the 
avalanche correlation function integrated over the field H (as 
in figure pi] ), while the other three are measurements of the 
correlation along the three axis, the six face diagonals, and 
the four body diagonals. Avalanches involving more than four 
spins show no noticeable anisotropy: the critical point ap- 
pears to have spherical symmetry. The same result is found 
in 2 dimensions. 

Figure [25| a shows the integrated avalanche correlation 
curves collapse in 3 dimensions for L = 320 and R > R c . 
The exponent v is obtained from such collapses by ex- 
trapolating to R — R c (figure [25^)) as was done for other 
collapses. The exponent j3/v can be obtained from these 
collapses too, but it is much better estimated from the 
magnetization discontinuity covered below. The value of 
0/v, listed in Table |l| alongside all the other exponents, 
is derived from the magnetization discontinuity collapses 
only. 



We have also looked for possible anisotropies in the 
integrated avalanche correlation function in 2 and 3 di- 
mensions. The anisotropic integrated avalanche corre- 
lation functions are measured along "generalized diago- 
nals": one along the three axis, the second along the six 
face diagonals, and the third along the four body diag- 
onals. We compare the integrated avalanche correlation 
function and the anisotropic integrated avalanche corre- 
lation functions to each other, and find no anisotropies 
in the correlation, as can be seen from figure [2?j . 

4- Spanning Avalanches 

The critical disorder R c was defined earlier as the dis- 
order R at which an infinite avalanche first appears in 
the system, in the thermodynamic limit, as the disor- 
der is lowered. At that point, the magnetization curve 
will show a discontinuity at the magnetization M C (R C ) 
and field H C (R C ). For each disorder R below the critical 
disorder, there is one infinite avalanche that occurs at 
a critical field H C (R) > H C (R C ) while above R c 

there are only finite avalanches. This is the behavior for 
an infinite size system. In a finite size system far below 
and above R c the above picture is still true, but close to 
the critical disorder, as we approach the transition, the 
avalanches get larger and larger, and we expect that one 
of them will be on the order of the system size and span 
the system from one "side" to another in at least one di- 
rection. This avalanche is not the infinite avalanche; it is 
only the largest avalanche that occurs close to the critical 
point. If the system was larger, this avalanche would be 
non-system spanning. Such an avalanche (which spans 
the system) we call a spanning avalanche. 

In our numerical simulation, we find that for finite sizes 
L, there are not one but many such avalanches in 4 and 
5 dimensions (and maybe 3), and that their number in- 
creases as the system size increases. Figures ^(a-c) show 
the number of spanning avalanches as a function of dis- 
order i?, for different sizes and dimensions. In 4 and 5 
dimensions, the spanning avalanche curves become more 
narrow as the system size is increased. Also, the peaks 
shift toward the critical value of the disorder (4.1 and 5.96 
respectively), and the number of spanning avalanches 
at R c increases. This suggests that in 4 and 5 dimen- 
sions, for L — > oo, there will be one infinite avalanche 
below R c , none above, and an infinite number of span- 
ning avalanches at the critical disorder R c . (These span- 
ning avalanches arc infinite avalanches for L — > oo.) In 
3 dimensions, the results are not conclusive, which can 
be noticed from figure |27|a, but also from the value of 
the spanning avalanche exponent 9 — 0.15 ± 0.15 defined 
below (a value of implies only one infinite or spanning 
avalanche at R c as L — > oo). 
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FIG. 27. (a) Number of spanning avalanches N in 3 
dimensions, occurring in the system between H = — oo to 
H — oo, as a function of the disorder R, for linear sizes L: 
20 (dot-dashed), 40 (long dashed), 80 (dashed), 160 (dotted), 
and 320 (solid). The critical disorder R c is at 2.16. The 
error bars for each curve tend to be smaller than the peak 
error bar for disorders above the peak and larger for disorders 
below the peak. They are not given here for clarity. Note 
that the number of avalanches increases only slightly as the 
size is increased, (b) Number of spanning avalanches in 
4 dimensions. The critical disorder is 4.1. (c) Number 
of spanning avalanches in 5 dimensions. The critical 
disorder is 5.96. Both in 4 and 5 dimensions, the peaks grow 
and shift towards R c as the size of the system is increased, (d) 
Collapse of the spanning avalanche curves in 4 dimensions for 
linear sizes L = 20, 40, and 80. The exponents are 6 = 0.32 
and v = 0.89, and the critical disorder is R c = 4.10. The 
collapse is done using r = (R c — R)/R. 

In percolation, a similar multiplicity of infinite clusters 
p7|p8| (as the system size is increased) is found for di- 
mensions above 6 which is the upper critical dimension 
(UCD). The UCD is the dimension at and above which 
the mean field exponents are valid. Below 6 dimensions, 
there is only one such infinite cluster. The existence of 
a diverging number of infinite clusters in percolation is 
associated with the breakdown of the hyperscaling rela- 
tion above 6 dimensions. Since a hyperscaling relation 
is a relation between critical exponents that includes the 
dimension d of the system, it is always only satisfied up 
to and including the upper critical dimension. In our 
system, the upper critical dimension is also 6, but we 
find spanning avalanches in dimensions even below that. 
In a comment by Maritan et al. p9| , it was suggested 
that our system should satisfy the hyperscaling relation: 
dv — = X/a which is also the one found in percola- 
tion |28). But since our system has spanning avalanches 
below the upper critical dimension, this hyperscaling re- 
lation breaks down below 6 dimensions. Due to the ex- 
istence of many spanning avalanches near R c , the new 
"violation of hyperscaling" relation for dimensions 3 and 
above becomes HJlol: 

(56) 

where 9 is the "breakdown of hyperscaling" or spanning 
avalanches exponent defined below. One can check that 
our exponents in 3, 4, and 5 dimensions and mean field 
satisfy this equation (see Tables || and III). 

For the simulation, we define a spanning avalanche to 
be an avalanche that spans the system in one direction. 
We average over all the directions to obtain better statis- 
tics. Depending on the size and dimension of the system 
and the distance from the critical disorder, the number of 
spanning avalanches for a particular value of disorder R is 
obtained by averaging over as few as 5 to as many as 2000 
different random field configurations. We define the ex- 
ponent 9 such that the number TV of spanning avalanches, 
at the critical disorder R c , increases with the linear sys- 
tem size as: N ~ L e (9 > 0). The finite size scaling 



(d-9)v-l3= X/a 
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form |l8|jr9| for the number of spanning avalanches close 
to the critical disorder is: 



N(L,R) ~ L e tf ± (LV v \r\) 



(57) 



where v is the correlation length exponent and Af± is the 
corresponding scaling function (± indicates the sign of 
r). The corrections to scaling are subdominant as ex- 
lained in appendix A. The collapse is shown in figure 
The values for 9 and v from collapses of curves of 
sizes L = 20,30,40, and 80 in 4 dimensions, are shown 
in Table IV. (We show the results and collapses in 4 di- 
mensions here since the existence of spanning avalanches 
in 3 dimensions is not conclusive.) These values are used 
along with the results from other collapses to obtain Ta- 
ble [n| In the analysis of the avalanche size distribution, 
magnetization, and correlation functions for R > i? c , how 
close we chose to come to the critical disorder R c was de- 
termined by the spanning avalanches: we include no val- 
ues R below the first value which exhibited a spanning 
avalanche. 



5. Magnetization Discontinuity 

We have mentioned earlier that in the thermodynamic 
limit, at and below the critical disorder i? c , there is a crit- 
ical field H C (R) > H C (R C ) at which the infinite avalanche 
occurs. Close to the critical transition, for r small and 
H = H C (R), the change in the magnetization due to the 
infinite avalanche scales as (equation (|3£ 



AM{R) 



(58) 



where r = (i? c — R)/R, while above the transition, there 
is no infinite avalanche. 




-rL 



FIG. 28. (a) Change in the magnetization due to the 
spanning avalanches in 4 dimensions, for several linear 
sizes L, as a function of the disorder R. (b) Scaling collapse 
of the curves in (a) using r = (R c — R)/R. The exponents 
are 1/u = 1.12 and j3/v = 0.19, and the critical disorder is 
R c = A.l. 

In finite size systems, the transition is not as sharp: we 
have spanning avalanches above the critical disorder. If 
we measure the change in the magnetization due to all the 
spanning avalanches (the infinite avalanche is included 
too), the scaling form for that quantity is going to de- 
pend on the system size L analogous to the scaling of the 
number of spanning avalanches: 



AM(L. R) ~ \rf AM±(L 1 ^\r\) 



(59) 



where AM± is a universal scaling function. (Since 
AM(L, R) is measured at h' = 0, corrections to scal- 
ing are subdominant; see also appendix A.) Defining a 
new universal scaling function AAA ± : 

AM±(L 1 ^\r\) = (L 1 / l/ |r|) _/3 AM±{L 1 ^ u \r\) (60) 

we obtain the scaling form: 

AM(L,R)~ L~PI U AM±{L 1 / V \r\) (61) 



Figures [28| a and |28p show the change in the magneti- 
zation due to the spanning avalanches in 4 dimensions, 
and a scaling collapse of that data (similar results exist 
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in 3 and 5 dimensions). Notice that as the system size 
increases, the curves approach the \r\@ behavior. The ex- 
ponents!/ v and (3/v are extracted from scaling collapses 
(figure p8p) and are listed in table 0. The value of (3 is 
calculated from (3/v and the knowledge of v, and is the 
value used for collapses of the magnetization curves (see 
earlier). 
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FIG. 29. (a) Second moments of the avalanche size dis- 
tribution integrated over the field H, in 5 dimensions. 

Error bars are largest for smaller disorders (shown on the 
curves). The curves have between 24 and 50 points, and 
the value of the second moment for each disorder is aver- 
aged over 3 to 100 different random field configurations, (b) 
Scaling collapse of the L = 10, 20, and 30 curves from (a) 
using r = (R c — R)/R. The exponents are \ jv — 1.47 and 
p — — (t + a (35 — 3)/av — 2.95, and the critical disorder is 
R c = 5.96. 



6. Moments of the Avalanche Size Distribution 

The second moment of the avalanche size distribution 
was defined earlier (see the mean field simulation sec- 
tion) . Wc found that the scaling form of the integrated 
over H second moment is (equation pq) : 



(S 2 ) int ~L-^-V/™S?\L^\r\) 



(62) 



where L is the linear size of the system, r is the reduced 
disorder, iS^_ is the scaling function, and v is the cor- 



relation length exponent. The corrections are subdomi- 
nant (appendix A). We can similarly define the third and 
fourth moment, with the exponent — (r + <j(35 — 3)/av 
replaced by — (t + a (35 — A)/<ji> and — (t + a (35 — b)/av 
respectively. Figures p9| a and p9| b show the second mo- 
ments data in 5 dimensions for sizes L = 5, 10, 20, and 
30, and a collapse (again, results in 3 and 4 dimensions 
are similar and we have chosen to show the curves in 
5 dimensions for variety). The curves are normalized 
by the average avalanche size integrated over all fields 
H: f*™f™S D{S,R,H,L) dS dH. The spanning 
avalanches and the infinite avalanche are not included 
in the calculation of the moments. The collapse does 
not include the L = 5 curve because, due to finite size ef- 
fects, this curve does not collapse well with the larger size 
curves. Tabic VI shows the values of the exponents and 
R c from the collapses. The exponents for the third and 
fourth moment can be calculated from this table, and we 
find that they agree with the values obtained from their 
respective collapses. 



7. Avalanche Time Measurement 

The exponents we have measured so far are static scal- 
ing exponents: they do not depend on the dynamics of 
the model. If we measure the time an avalanche takes 
to occur, we are making a dynamical measurement. The 
time measurement in the numerical simulation is done 
by increasing the time "meter" by one for each shell of 
spins in the avalanche; it corresponds to a synchronous 
dynamics, where, when all unstable spins are flipped, 
time is incremented by one, and the new list of unsta- 
ble spins is generated. The scaling relation between the 
time t it takes an avalanche to occur and the size S of 
that avalanche for small disorder r can be found by not- 
ing that the characteristic duration of an avalanche is 
proportional to the correlation length £ to the power z 
MM: 



(63) 



The exponent z is known as the dynamical critical ex- 
ponent. Equation ( |63| ) gives the scaling for the time it 
takes for a spin to "feel" the effect of another a distance £ 
away. Since the correlation length £ scales like r~ v close 
to the critical disorder, and the characteristic size S as 
r -1 /' 7 , the time t then scales with large sizes as: 



t 



(64) 
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FIG. 30. (a) Avalanche time distribution curves in 
3 dimensions, for avalanche size bins from about 2000 to 
40000 spins (from upper left to lower right corner). The sys- 
tem size is 800 3 at R = 2.26. The curves are from only one 
random field configuration, (b) Scaling collapse of curves 
in (a). The values of the exponents are avz = 0.57 and 
(r + a/36 + avz) I avz = 4.0. 

In our simulation, we measure the distribution of times 
for each avalanche size S. The distribution of times 
D t (S, R, H, t) for an avalanche of size S close to the crit- 
ical field H c and critical disorder R c is 



D t (S,R, H,t) ~ S-i V ( i ) (S a \r\,h/\r\ 0S ,t/S' J " 
where q = r + avz, and is defined such that 

/+oo />oo 
/ D t {S,R,H,t) dH dt = 
-oo Jl 

S -(T + *0S) p(™*)(^| r |) 



(65) 



(66) 



where p^ n *^ was defined in the integrated avalanche size 
distribution section. The avalanche time distribution in- 
tegrated over the field H, at the critical disorder (r = 0) 
is: 



D[ mt) (S,t) 



j — (T+al3&+avz)/(jvz jy(int) Ujgwz ^ (67) 



which is obtained from equation ( pq ) in a derivation anal- 
ogous to the one for the integrated avalanche size distri- 
bution scaling form. 



Figures |30| a and pOp show the avalanche time distri- 
bution integrated over the field H for different avalanche 
sizes, and a collapse of these curves using the above scal- 
ing form, for a 800 3 system at R = 2.260 (just above 
the range where spanning avalanches occur). The data 
is saved in logarithmic size bins, each about 1.2 times 
larger than the previous one. The time is also measured 
logarithmically (next bin is 1.1 times larger than the pre- 
vious one). The extracted value for z in 3 dimensions is 
1.68 ± 0.07. The results for other dimensions are listed 
in Table [f| 

C. Simulation Results in 2 Dimensions 

The critical transition in the shape of the hysteresis 
loop is observed in the simulation, and expected from the 
renormalization group |l^]l5| |, in 3, 4, and 5 dimensions. 
We also found that the upper critical dimension, at and 
above which the mean field exponents become correct, is 
six. Furthermore, in one dimension, we expect that in a 
thermodynamic system, with an unbounded distribution 
of random fields, there will be no infinite avalanche for 
R > 0. That will be so because if there is any random- 
ness, there will be a spin in the linear chain that will 
have the "right" value for its random field to stop the 
first avalanche. For a bounded distribution of random 
fields, the scaling behavior near the transition will not 
be universal H|; instead, it will depend on the exact 
shape of the tails of the distribution of random fields. 
Then, the question that remains is: what happens in two 
dimensions? 

From the simulation and a few arguments that we are 
about to show, we conjecture that the two dimensional 
exponents will have the values: r + crj38 = 2, r = 3/2, 
1/^ = 0, and av = 1/2. (The other exponents (except z) 
can be found from exponent relations [l4|Jl^| using these 
values.) The "arguments" are as follows. 

It is quite possible that two is the lower critical dimen- 
sion (LCD) for our system. At the lower critical dimen- 
sion, the critical exponents are often ratios of small in- 
tegers, and it is often possible to derive exact solutions. 
Since the geometry in 2 dimensions allows for at most 
one system spanning avalanche, the "breakdown of hy- 
perscaling" exponent 9 (see section IV B) must be zero, 
and the hyperscaling relation HJlq] is restored: 



(68) 



av v 

We know that this relation is violated in 4 and 5 di- 
mensions, and is probably violated in 3 dimensions. In 
dimensions above two, the hyperscaling relation is mod- 
ified by the exponent 9 which gives a measure of the 
number of spanning avalanches near the critical transi- 
tion, as a function of the system size. Figure |3l] shows 
the number of spanning avalanches in 2 dimensions for 
several system sizes. Notice that, as assumed, there is 
not more than one spanning avalanche in the system. 
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FIG. 31. Number of spanning avalanches in 2 dimen- 
sions as a function of disorder R, for several system sizes. The 
data points are averages between as little as 10 to as many 
as 2200 random field configurations. Some typical error bars 
near the center of the curves are shown; error bars are smaller 
toward the ends. Note that there is no more than one span- 
ning avalanche. 

We use two more arguments to derive the critical ex- 
ponents. In 2 dimensions, we find that the avalanches 
"look" compact (figure |32|). (The avalanches in 3 dimen- 
sions are not compact (figure [53]).) This implies that 
1/av = d = 2, which leads to (3/v = from equation 
(pq). Furthermore, it is often the case that in the lower 
critical dimension, the Harris criterion 



v 2 
— > - 

" d 



35 



(69) 



becomes saturated (an equality); so in 2 dimensions we 
expect (36 /u = 1. From this and the previous result, the 
exponent which gives the decay in space of the avalanche 
correlation function 



9 ^/3 (36 

r] = 2 H 

v v 



(70) 



(see references for the derivation of all the expo- 

nent relations) becomes equal to r\ = 1. 

Since at the LCD the correlation length typically di- 
verges exponentially as the critical point is approached, 
we expect v — > oo, and (3 can be finite. Using the expo- 
nent relation @^6|: 

T-2 = a/3(l-6), (71) 

we further find that r = 3/2 and r + a (3d = 2. 




FIG. 32. Simulation in 2 dimensions of a 400 spin system 
at R — 0.8. The figure shows the configuration of the system 
after a spanning avalanche has just occurred (grey region). 
The dark area corresponds to spins that have not yet flipped, 
while the white area are spins that have flipped earlier. Notice 
that the spanning avalanche (grey area) seems compact. 

We must mention that our firm conjectures about the 
exponents in two dimensions must be contrasted with 
our lack of knowledge about the proper scaling forms. 
As mentioned above, at the LCD the correlation expo- 
nent v typically diverges, although some combinations of 
critical exponents stay finite (hence av = 1/2). Those 
which diverge and those which go to zero usually must 
be replaced by exponents and logs, respectively. We have 
used three different RG-scaling ansdtze to model the data 
in two dimensions. (1) We used the traditional scal- 
ing form £ ~ \R C — R\~ v , deriving v — 5.3 ± 1.4 and 
R c = 0.54 ± 0.04. These collapses worked as well as 
any, but the large value for v (and larger value still for 
1/(7 = 10±2) makes one suspicious. (2) We used a scal- 
ing form suggested by Bray and Moore |3^] in the con- 
text of the equilibrium thermal random field Ising model 
at the LCD, where R c = 0: if they assume that R is 
a marginal direction, then by symmetry the flows must 
start with i? 3 , leading to £ - e (<V|Re--R| 2 ) = e (a/R 2 )^ 
This form has the fewest free parameters, and most of 
the collapses were about as good as the others (except 
notably for the finite-size scaling of the moments of the 
avalanche size distribution, which did not collapse well 
once spanning avalanches became common). (3) We de- 
veloped another possible scaling form, based on a finite 
R c and R marginal, which generically has a quadratic 
flow under coarse-graining: here £ ~ e { b /\R<=-R\) _ ^ye 
find R c = 0.42 ± 0.04. The rational behind these three 
forms is shown in appendix A. 

The results from data collapses in two dimensions were 
obtained from measurements of the spanning avalanches, 
the second moments of the avalanche size distribu- 
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tion, the integrated avalanche size distributions, and the 
avalanche correlations. The magnetization curves are 
also obtained from the simulation, but as in the higher 
dimensions, the scaling region is small (around H c and 
M c ), and the collapses do not define the exponents well. 



using the three different scaling forms are shown in figures 
|J(b-d). The first one (figure fjjb) is: 




FIG. 33. Largest avalanche occurring in the hysteresis loop 
in a 40 3 spins system near the critical point. The avalanche 
is not compact. 

Measurements that require the knowledge of the criti- 
cal randomness are the binned avalanche size distribution 
from which we extract the exponents r and f38, the crit- 
ical magnetic field iJ c , and the avalanche time measure- 
ment which gives the exponent z. These measurements 
were not obtained at the critical disorder because R c is 
not well defined as was mentioned above, and because 
for low disorders (less than 0.71 for a 7000 2 system), the 
system flips in one infinite avalanche, and such measure- 
ments are therefore not possible. We have nevertheless 
estimated the values of some of these exponents and of 
H c , from data obtained at a larger disorder (where there 
is no spanning avalanche). From the avalanche size dis- 
tribution binned in H at R — 0.71 and L = 7000, and 
the magnetization curves, we find that the critical field 
H c is around 1.32. A straight line fit through the data 
agrees with a possible value of r = 3/2 (the conjectured 
value). From the time distribution of avalanche sizes for 
a system of 30000 2 spins, at R — 0.65, we measured (from 
a straight linear fit) the exponent avz to be 0.64. The 
other exponents were obtained from scaling collapses as 
follows. 

Figure |34| a shows the second moments of the avalanche 
size distribution for several system sizes. The collapses 



(S 2 ) int ~L-^ s - 3 »™ S™( L \r\") 



(72) 



which is the kind of scaling form used in 3, 4, and 5 
dimensions. This form assumes £ ~ |r| ~ v . The expo- 
nents are (r + af38 — 3)/av — —1.9 and v = 5.25, and 
r = (R c — R)/R with R c = 0.54. The second scaling form 
(figure [54|c) is: 

(S 2 ) mt ~ L-^'-»/- S%(L e -a/l*-*l 2 ) (73) 

which is obtained from £ ~ e O/l-Rc--R| )_ The vauies f 
the exponents and parameters are: (r + a (38 — Z)/av = 
— 1.9, a — 3.4 (5 is not universal), and R c = (by as- 
sumption; see previous paragraph). Notice that this col- 
lapse is not as good as the other two; a better collapse is 
obtained with R = 0.15 and a = 2.0. If this is the correct 
scaling form and R c = 0, this discrepancy can be due to 
finite size effects. The third scaling form is (figure f34]d): 

(S 2 ) mt - L-C^- 3 )/- Sg>(L e-VI^-HI) (74) 

which is obtained from £ ~ e ( b /\ R c~ R \) . The values of the 
exponents and parameters are: (t + a(38 — 3)/ai> = —1.9, 
b = 2.05 (b is also non-universal), and R c — 0.42. As it 
is clear from the last three figures, collapses with these 
different scaling forms are comparable. Notice that the 
exponent (r + a j38 — i)/av is the same for the three 
collapses, but that \jv is zero for the last two (by as- 
sumption) while it is 0.19 for the first collapse. Let's 
now look at the collapses of the integrated avalanche size 
distribution curves, which are not finite size scaling mea- 
surements. 
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FIG. 34. (a) Second moments of the avalanche size 
distribution in 2 dimensions, integrated over the external 
field H, for several system sizes. The data points are aver- 
ages over up to 2200 random field configurations. Error bars 
are smaller than shown for larger disorders, (b), (c), and (d) 
Scaling collapses of the second moments of the avalanche size 
distribution in 2 dimensions, integrated over the field H. The 
curves that are collapsed are of size: 50 2 , 100 2 , 300 2 , 500 2 , 
1000 2 , 3000 2 , 5000 2 , 7000 2 , and 30000 2 . See text for the scal- 
ing forms, and the values of the exponents and parameters. 



Figure |35| a shows the integrated avalanche size distri- 
bution curves for a 7000 2 spin system, at several values 
of the disorder R. Earlier, in figure we saw the fit to 
the scaling collapse of such curves, done using the same 
scaling form as in 3, 4, and 5 dimensions: 



D int (S,R) ~ S-^+'W V ( " lt) {S a \r\) 



(75) 



(The — sign indicates that the collapsed curves are for 
r < 0, ic. R> R c .) However, S a \r\ might not be the ap- 
propriate scaling argument in 2 dimensions. First, from 
figure |2(], the scaling curve in 2 dimensions differs dra- 
matically from the scaling curves in higher dimensions 
for small arguments X — S a \r\. The mean field scal- 
ing function -D ( ' nt) {X) is a polynomial for small X, and 
we expected (and found) a similar behavior in 5, 4 and 
3 dimensions (but notice that the scaling function in 3 
dimensions is starting to look like the curve in 2 dimen- 
sions for small X). In 2 dimensions, if we collapse our 
data (figure |35|b) using the scaling form: 



D mt (S,R) - S-^ +a ^ V ( -i nt \s\r\ 1 ' a ) (76) 

with t + a (35 = 2.04, 1/a = 10, and r = (R c - R)/R, 
we find that the scaling function for small X = S'|r| 1 / cr 
looks linear with power one! This might imply that the 
scaling function D 1 ^ 71 ^ (S'|r| 1 /' J ) (eqn. (|76|)) is the one that 
is analytic for small arguments in 2 dimensions. 
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FIG. 35. (a) Integrated avalanche size distribution 

curves for several disorders in 2 dimensions, at the system 
size L = 7000. The curves are averages over 10 to 20 ran- 
dom field configurations, and have been smoothed, (b), (c), 
and (d) Scaling collapses of the data from (a) using the three 
scaling forms and the exponents from the text. The collapsed 
curves have disorders: 0.72, 0.74, 0.77, and 0.80. The straight 
grey line in each of the plots has a slope of one. 

Second, we conjectured above that the values for 
a and \jv are probably zero in 2 dimensions, 
and that only the combination o~v is finite (av — 
1/2). It follows that, for the other two scaling 
forms we use, the arguments of the scaling function 
should be and Se -b/<rv\R-R c \ ^ and not 

S a e -a/v\R-R a \ 2 and ga e -b/v\R-R c \ respectively. This is 
analogous to using 5'|r| 1 / CT in the scaling form (|7^). We 
should mention here that both equation (]75|) and equa- 
tion ( f76|) give the same scaling exponents T + af36 and a, 
and that in all our scaling collapses we have assumed that 
the same scaling argument is valid for small and large X 
(and in between). This in general, does not have to be 
true. 




Avalanche Size (S) 
FIG. 36. Integrated avalanche size distribution 
curve in 2 dimensions, for a system of 30000 2 spins, at 
R — 0.650. Shown are two linear fits to the data: one for 
small sizes and the other for large sizes. The slope for the 
fit at small S is 0.90. The fit was done for sizes in the range 
[10, 250]. The slope differs by less than 5% when the range is 
changed (5* is never larger than 400 though) . The slope for 
the fit at large 5* is 1.78. The slope differs by less than 2% 
when the range is changed (S is never smaller than 10000). 
The conjectured value for r + a (38 is 2 which is different from 
1.78. This is similar to the behavior we saw in 3, 4, and 5 
dimensions. On the other hand, for small sizes we expect the 
exponent T + a/3S— 1 — 1 (see text). Again, the two measure- 
ments don't completely agree, but the slope from our data 
does seem to indicate such a behavior. 

Equation ( f76| ) is therefore one of the three scaling forms 
we use. The second scaling form is: 
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D int (S,R) ~ S 



HI 1 



(77) 



shown in figure p5p , with r+afid = 2.04, a/av = 7.0 (this 
implies that av = 0.49), and r — R c — R with R c = by 
assumption. And finally, the third scaling form we use 



D int (S,R) ~ ^(int)(l) / 5e -6/^|fl c -flA 



(78) 



shown in figure |35|d, with r + a/36 = 2.04, &/W = 4.0 
(which makes av — 0.51), and r = R c —R with R c = 0.42. 
Again, not only are all three collapses comparable, but 
the exponents extracted from them are as well. The ex- 
ponent for the slope of the distribution is r + a (36 = 2.04 
for the three collapses, and the exponent combination 
av is around 0.51 (for the first collapse a = 0.10, while 
v = 5.25 from the equivalent second moment collapse). 

Figures ^(b-d) show that the scaling function 2?^™'' 
seems to be linear with slope one for small arguments 
(the grey lines have slope one) and that the constant 
term in the polynomial expansion is zero (or close to 
zero) . This leads to a singular scaling function correction 
to the avalanche size distribution exponent r + a[35 for 
small non-zero X: 

(79) 

(Note that we could have used £>f m *-"- 1 ' ) or p^"*^ 2 ) as 
well.) 

Recall that because of the "bump" in the avalanche 
size distribution scaling function in 3, 4, and 5 dimen- 
sions, and in mean field, the slope of the raw data curves 
did not agree with the value of the exponent t + a(35. 
In 2 dimensions, this is still true, but we also find a sin- 
gular behavior for D^ nt \x), which changes the slope of 
the data curve for small X . In figure |36|, an integrated 
avalanche size distribution curve for a system of 30000 2 
spins, at R = 0.65, is plotted along with the linear fits 
to the data for small and large size S. For large S, the 
slope is close to but not equal to 2, while for small S, the 
slope is close to one! 

The avalanche correlation data (see figure |37|a) is col- 
lapsed with three different scaling forms as well. These 
forms are analogous to the ones used for the second mo- 
ments collapses, but with the distance x taking the place 
of the system size L. The collapses and the extracted 
exponents from these three forms are again very simi- 
lar, and only one of the collapses is shown in figure |37|b. 
The value of j3/v from these collapses is 0.03 ± 0.06. If 
we compare figure |37j b with the collapse of the avalanche 
correlation in 3 dimensions (fig. |25|a), we find that the 
scaling function in 2 dimensions seems to be singular 



with slope one for small distances, as is the integrated 
avalanche size distribution for small sizes. 

The spanning avalanches data are also analyzed using 
three scaling forms similar to those used for the second 
moments of the avalanche size distribution collapses. The 
exponent 9 is poorly defined from these collapses (and is 
therefore not listed in Table VII), although the data does 
collapse for the exact value: 8 = 0. 

The three collapses for all the measurements we have 
done are very similar. This is not a surprise: it is always 
hard to distinguish large power laws (v and l/a are large 
in the "linear argument" scaling form (eqns. ( ]7^ ) and 
([76]))) from exponentials. Although some of the expo- 
nents have very different values in the three collapses, the 
average of the exponents from the three methods agree 
within the error bars with each method (see figure |38| ) 
and our conjectures. In conclusion, although we do not 
know the correct scaling form for the data in 2 dimen- 
sions, the possible three scaling forms we mention give 
exponent values that are compat ible with eac h othe r and 
with our conjectures (see Table VII). (Table VIII gives 
the conjectured values for the exponents that have not 
been measured in the collapses.) Much larger system 
sizes might be necessary to obtain more conclusive re- 
sults. 
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FIG. 37. (a) Avalanche correlation function in 2 di- 
mensions, integrated over the external field H , for several 
disorders R and the system size L — 30000. Only the curve 
with the smallest disorder is an average over several random 
field configuration, (b) Scaling collapse of the avalanche cor- 
relation curves in 2 dimensions, for a system of 30000 2 spins. 
The exponent values are: v = 5.25 and j3/v = 0. The critical 
disorder is R c = 0.54, and r = (R — R c )/R. Notice that for 
small aj|r| 1 ', the scaling function looks singular with a power 
close to one (the straight line has a slope of one). 



V. COMPARISON WITH THE ANALYTICAL 
RESULTS 

We have compared the simulation results with the 
renormalization group analysis of the same system 
p4] , p^| . According to the renormalization group the up- 
per critical dimension (UCD), at and above which the 
critical exponents are equal to the mean field values, is 
six. Close to the UCD, it is possible to do a 6 — e expan- 
sion (e is small and greater than 0) , and obtain estimates 
for the critical exponents and the magnetization scaling 
function, which can then be compared with our numerical 
results. Furthermore, at dimension eight there is a pre- 
diction for another transition. Below eight dimensions, 
there is a discontinuity in the slope of the magnetization 
curve as it approaches the "jump" in the magnetization 
(R < R c ), while above eight dimensions the approach is 
smooth. 
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FIG. 38. Numerical values (filled symbols) of the expo- 
nents t + af35, r, 1/v, avz, and av (circles, diamond, tri- 
angles up, squares, and triangle left) in 2, 3, 4, and 5 dimen- 
sions. The empty symbols are values for these exponents in 
mean field (dimension 6). Note that the value of r in 2d is 
the conjectured value: we have not extracted r from scaling 
collapses (see text). We have simulated sizes up to 30000 2 , 
1000 3 , 80 4 , and 50 5 , where for 320 3 for example, more than 
700 different random field configurations were measured. The 
long-dashed lines are the e expansions to first order for the 
exponents r+a/35, r, avz, and av. The short-dashed lines are 
Borel sums |j3| for \ jv. The lowest is the variable-pole Borel 
sum from LeGuillou et al. the middle uses the method 
of Vladimirov et al. to fifth order, and the upper uses the 
method of LeGuillou et al., but without the pole and with 
the correct fifth order term. The error bars denote system- 
atic errors in finding the exponents from extrapolation of the 
values obtained from collapses of curves at different disorders 
R. Statistical errors are smaller. 

Figure ^ shows the numerical and analytical results 
for five of the critical exponents obtained in dimensions 
two to six (in six dimensions, the values are the mean field 
ones). The other exponents can be obtained from scaling 
relations uA 



3N 



arc 



3 4 5 
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Jl6|| . The exponent values in figure 

obtained by extrapolating the results of scaling collapses 
to either R ^ R c or l/i^0 (see section on simulation 
results). In two dimensions, which is possibly the lower 
critical dimension, the plotted values are averages from 
the three different scaling forms used to collapse the data 
and extract the exponents. The error bars shown span all 
three ansatze, and are compatible with our conjectures 
from the previous section. The long-dashed lines are the e 
expansions to first order for T + a(3S, r, avz, and av. The 
three short-dashed lines Jl4|] are Borel sums |33| for \ jv. 
The lowest is the variable-pole Borel sum from LeGuillou 
et al. ]33|], the middle uses the method of Vladimirov 
et al. to fifth order, and the upper uses the method 
of LeGuillou et al., but without the pole and with the 
correct fifth order term |14|. Notice that the numerical 
values converge nicely to the mean field predictions, as 
the dimension approaches six, and that the agreement 
between the numerical values and the e expansion is quite 
impressive. 
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FIG. 39. Comparison between six simulation curves (thin 
lines) and the dM/dH curve (thick dashed line) obtained from 
a parametric form |}5| to third order in e. The six curves 
are for a system of 30 5 spins at disorders: 7.0, 7.3, and 7.5 
(R c = 5.96 in 5 dimensions), and for a system of 50 5 spins at 
disorders: 6.3, 6.4, and 6.5 (for larger fields, these are closer 
to the dashed line in the figure). All the curves have been 
stretched/shrunk in the horizontal and vertical direction to 
lie on each other, and shifted horizontally. 

The e expansion can be an even more powerful tool if 
it can predict the scaling functions. This has been done 
for the magnetization scaling function of the pure Ising 
model in 4 — e dimensions |34|,|35| . Since the e expansion 
for our model is the same as the one for the equilibrium 
RFIM (TJ] , and the latter has been mapped to all orders 
in e to the corresponding expansion of the regular Ising 
model in two lower dimensions |l4|,|36|,[3^] , we can use the 
results obtained in J34p5| ]. This is done in figure |3^, 
which shows the comparison between the dM/dH curves 
obtained in 5 dimensions at R = 6.3, 6.4, 6.5, 7.0, 7.3, 7.5 
(R c = 5.96) (the curves have been stretched/shrunk to 
lie on top of each other, and shifted horizontally so that 
the peaks align), and the parametric form (thick dashed 
line) for the scaling function of dM/dH, to third order 
in e, where e = 1 in 5 dimensions (see pSfl). As we see, 
the agreement is very good in the scaling region (close to 
the peaks). This brings up the possibility of using the e 
expansion for the scaling function to extract the critical 
exponents from simulation or experimental data. So far 
though, only the scaling function for the magnetization 
has been obtained. 



-1.00 ' ■ 1 ■ 1 ■ ' 

1.0 3.0 5.0 7.0 

Magnetic Field (H/J) 
FIG. 40. Magnetization curves in 5, 7, and 9 dimen- 
sions. The disorders for these curves are R = 3.3, 4.7, and 
6.0 for 30 5 , 10 7 , and 5 9 size systems respectively. The dashed 
lines represent the "jump" in the magnetization. Notice that 
in 9 dimensions the approach to the "jump" seems to be con- 
tinuous. 

As another check between the simulation and the 
renormalization group, we have looked for the predicted 
transition in eight dimensions. Figure shows the mag- 
netization curves in 5, 7, and 9 dimensions (system sizes: 
30 5 , 10 7 , and 5 9 ) for values of the disorder equal to |d, 
where d is the dimension. These values of disorder are 
below the critical disorder in dimensions below six, and 
are expected to be below for dimensions 7 and 9 as well. 
For 5d and 7d, the approach to the "jump" in the mag- 
netization is discontinuous. Above the eight dimension, 
the approach is continuous (see close ups in figure fill). 
This is as expected from the renormalization group [ [14] . 
We have also looked at dM/dH. which appears clearly 
to diverge in d = 9 and not in d = 7 (figure |42| ) . 

VI. CONCLUSION 

We have used the zero temperature random field Ising 
model, with a Gaussian distribution of random fields, 
to model a random system that exhibits hysteresis. We 
found that the model has a transition in the shape of the 
hysteresis loop, and that the transition is critical. The 
tunable parameters are the amount of disorder R and 
the external magnetic field H. The transition is marked 
by the appearance of an infinite avalanche in the ther- 
modynamic system. Near the critical point, (Rc, He), 
the scaling region is quite large: the system can exhibit 
power law behavior for several decades, and still not be 
near the critical transition. This is important to keep in 
mind whenever experimental data are analyzed. If a tun- 
able parameter can be found, a system that appears to 
be SOC, might in reality have a disorder induced critical 
point. 
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FIG. 41. (a) and (b) Closeup of the magnetization 
curves in 7 and 9 dimensions respectively from figure 
In 8 dimensions, there is a prediction from the renormaliza- 
tion group jll|] that there is a transition in the way the jump 
is approached (see text). 

We have extacted critical exponents for the magneti- 
zation, the avalanche size distribution (integrated over 
the field and binned in the field), the moments of the 
avalanche size distribution, the avalanche correlation, the 
number of spanning avalanches, and the distribution of 
times for different ava lanc he sizes. These values are listed 
in Table [n] and Table VII, and were obtained as an aver- 
age of the extrapolation results (to R — > R c or L — > oo) 
from several measurements. For example, the correlation 
length exponent v is the average value from three dif- 
ferent collapses: the correlation function, the spanning 
avalanches, and the second moments of the avalanche 
size distribution, while the critical disorder R c is esti- 
mated from both the spanning avalanches collapses and 
the collapses of the moments of the avalanche size dis- 
tribution. As shown earlier, the numerical results com- 
pare well with the e expansion MJTq]. Furthermore, the 



renormalization group work predicts another transition 
in eight dimensions, which we find in the simulation as 
well. Comparisons to experimental Barkhausen noise 
measurements [ [l2| are very encouraging, and a more 
comprehensive review of possible experiments that ex- 
hibit disorder-driven critical phenomena similar to our 
model is under way p6[ . 
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FIG. 42. Derivative of the magnetization with respect to 
the field H, for the curves in figure [Io| The approach to the 
"infinite jump" seems to be continuous in 9 dimensions. Note 
that the vertical axis is logarithmic. 

Finally, we should mention that there are other mod- 
els for avalanches in disordered magnets. There is a large 
body of work on depinning transitions and the motion of 
the single interface [ p8|j3S| ] . In these models, avalanches 
occur only at the growing interface. Our model though, 
deals with many interacting interfaces: avalanches can 
grow anywhere in the system. Similar models exist with 
random bonds |^,^l| and random anisotropics. In the 
random bonds model, the interaction between neigh- 
boring spins i and j is random. 

The zero temperature random bond Ising model [}l0|j4l[ 
also exhibits a critical transition in the shape of the hys- 
teresis loop, where the mean bond strength is analogous 
to our disorder R. It has been argued numerically Glj] 
and analytically fljij , that as long as there are no long- 
range forces H and correlated disorder, the random bond 
and the random field Ising model are in the same uni- 
versality class. A comparison between our simulation 
and the results in reference 0j show that the 3 dimen- 
sional results agree quite nicely. However, in 2 dimen- 
sions, there are large differences, which we believe occur 
because of the small system sizes used by the authors for 
their simulation (only up to L — 100). We have seen that 
our results (see section on the 2 dimensional simulation) 
are very size dependent. Looking back for example at 
figure |3l], we find that for a system of L = 100 spins, a 
"good" estimate for the critical disorder R c would indeed 
be 0.75 as was found in @. However, we find after in- 
creasing the system size that the critical disorder R c is 
0.54 or lower. 
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APPENDIX A: DERIVATION OF THE VARIOUS 
SCALING FORMS AND CORRECTIONS 

In this paper we make extensive use of scaling col- 
lapses. Many variations are important to us: Widom 
scaling, finite-size scaling, singular corrections to scaling, 
analytic corrections to scaling, rotating axes, and expo- 
nentially diverging correlation length scaling. The un- 
derlying theoretical framework for scaling is given by the 
renormalization group, developed by Wilson and Fisher 
|j42f in the context of equilibrium critical phenomena and 
by now well explicated in a variety of texts Jl8|,|3"l 43 1. 

We have discovered that we can derive all the scaling 
forms and corrections that have been important to us 
from two simple hypotheses (found in critical regions): 
universality and invariance under reparameterizations. 
Universality is the statement that two completely dif- 
ferent systems will behave the same near their critical 
point [H] (for example, they can have exactly the same 
kinds of correlations). Reparameterization invariance is 
the statement that smooth changes in the units or meth- 
ods of measurement should not affect the critical proper- 
ties. We use these properties to develop the scaling forms 
and corrections we use in this paper. Each example we 
cover will build on the previous ones while developing a 
new idea. 

For our first example, consider some property F of a 
system at its critical point, as a function of a scale x. F 
might be the spin-spin correlation function as a function 
of distance x (or it might be the avalanche probability 
distribution as a function of size x, etc.) If two different 
experimental systems are at the same critical point, their 
.F's must agree. It would seem clear that they cannot be 
expected to be equal to one another: the overall scale 
of F and the scale of x will depend on the microscopic 
structure of the materials. The best one could imagine 
would be that 



F^xi) =AF 2 (Bx 2 ) 



(Al) 



where A would give the ratio of, say, the squared mag- 
netic moment per domain of the two materials, and B 



gives the ratio of the domain sizes. 

Now, consider comparing a system with itself, but with 
a different measuring apparatus. Universality in this self- 
referential sense would imply F{x) — AF(Bx), for suit- 
able A and B. If instead of using finite constant A and 
B, we arrange for an infinitesimal change in the measure- 
ment of length scales, we find: 



F(x) = (1-cee) F((l-e)x^ 



(A2) 



where e is small and a is some constant. Taking the 
derivative of both sides with respect to e and evaluating 
it at e = 0, we find — aF = xF' , so 



F(x) 



(A3) 



The function F is a power-law! The underlying reason 
why power-laws are seen at critical points is that power 
laws look the same at different scales. 

Now consider a new measurement with a distorted 

where 



F[B(x) 



measuring apparatus. Now F{x) ~ A 

A and B are some nonlinear functions. For example, 
one might measure the number of microscopic domains x 
flipped in an avalanche, or one might measure the total 
acoustic power B(x) emitted during the avalanche; these 
two "sizes" should roughly scale with one another, but 
nonlinear amplifications will occur while the spatial ex- 
tent of the avalanche is small compared to the wavelength 
of sound emitted: we expand B(x) = Bx + bo + bi/x+ . . . 
Similarly, our microphone may be nonlinear at large 
sound amplitudes, or the absorption of sound in the 
medium may be nonlinear: A(F) = AF + a 2 F 2 + . . . 
So, 



A 



F 



(B(x)) 



A^F(Bx)+F'{Bx)(b + b 1 /x- 



F"(Bx)(. 



a 2 F 2 (Bx) 



(A4) 



We can certainly see that our assumption of universality 
cannot hold everywhere: for large F or s mall x the as- 
sumption of reparameterization invariance ( |A4|) prevents 
any simple universal form. Where is universality pos- 
sible? We can take the power-law form F(x) ~ x~ a = 
x iog A/ log b whjch j s th e only for m all owed by linear repa- 
rameterizations and plug it into (A4), and we see that all 
these nonlinear corrections are subdominant (i.e., small) 
for large x and small F (presuming a > 0). If a > 1, 
the leading correction is due to bo and we expect 
corrections to the universal power law at small distances; 
if < a < 1 the dominant correction is due to a 2 , and 
we expect corrections of order x~ 2a . Thus our assump- 
tions of universality and reparameterization invariance 
both lead us to the power-law scaling forms and inform 
us as to some expected deviations from these forms. No- 
tice that the simple rescaling led to the universal power- 
law predictions, and that the more complicated nonlinear 
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rescalings taught us about the dominant corrections: this 
will keep happening with our other examples. 

For our second example, let us consider a property K 
of a system, as a function of some external parameter 
R, as we vary R through the critical point R c for the 
material (so r = R — R c is small). K might represent 
the second moment of the avalanche size distribution, 
where R would represent the value of the randomness; 
alternatively K might represent the fractional change in 
magnetization AM at the infinite avalanche ... If two 
different experimental systems are both near their critical 
points (ri and r 2 both small), then universality demands 
that the dependence of K\ and K 2 on "temperature" R 
must agree, up to overall changes in scale. Thus, using a 

simple linear rescaling K(r) = (1 — fiejK^Jl — e)rj leads 

as above to the prediction 



K{r) = r 



(A5) 



Now let us consider nonlinear rescalings, somewhat dif- 
ferent than the one discussed above. In particular, the 
nonlinearity of our measurement of K can be dependent 
on r. So, 



A 



(K{T) 



ao + a-iT + a 2 r 2 + . . . + a 01 K(r) + 



(A6) 



If fi > 0, these analytic corrections don't change the dom- 
inant power law near r = 0. However, if < 0, all the 
terms a n for n < —ji will be more important than the sin- 
gular term! Only after fitting them to the data and sub- 
tracting them will the residual singularity be measurable. 
For the fractional change in magnetization: AM ~ r@ 
has < < 1 (at least above three dimensions), so we 
might think we need to subtract off a constant term ao, 
but AM = for R > R c , so ao is zero. On the other 
hand, in a previous paper |Q, we discussed the singu- 
larity in the area of the hysteresis loop: Area ~ r 2_Q , 
where 2 — a = j3 + (35 is an analogue to the specific heat 
in thermal systems. Since a is near zero (slightly positive 
from our estimates of (3 and 5 in 3, 4, and 5 dimensions), 
measuring it would necessitate our fitting and subtract- 
ing three terms (constant, linear, and quadratic in r): we 
did not measure the area for that reason. 

For our third example, let's consider a function F(x, r), 
depending on both a scale x and an external parameter 
r. For example, F might be the probability Di nt that 
an avalanche of size x will occur during a hysteresis loop 
at disorder r = R — R c . Universality implies that two 
different systems must have the same F up to changes in 
scale, and therefore that F measured at one r must have 
the same form as if measured at a different r. To start 
with, we consider a simple linear rescaling: 



F 



(x,r) = (l-ae) f((1 - e)x, (1 + Ce)r). 



(A7) 



Taking the derivative of both sides with respect to e gives 
a partial differential equation that can be manipulated to 



show F has a scaling form. Instead, we change variables 
to a new variable y = x^r (which satisfies y' — y to order 
e). If F(x,y) = F(x,r) is our function measured in the 
new variables, then 

F(x, r) = F(x, y) = (1 - ae) f((1 - e)x, yj (A8) 

and —aF = xdF/dx shows that at fixed y, F ~ x~ a , 
with a coefficient T(y) which can depend on y. Hence 
we get the scaling form 



F(x,r) - x~ a T{x c r). 



(A9) 



This is just Widom scaling. The critical exponents a and 
C and the scaling function T{x^r) are universal (two 
different systems near their critical point will have the 
same critical exponents and scaling functions). We don't 
need to discuss corrections to scaling for this case, as 
they are similar to those discussed above and below (and 
because none were dominant in our cases). 

Notice that if we sit at the cri tica l point r — 0, the 
above result reduces to equation (A3) so long as .F(O) is 
not zero or infinity. If, on the other hand, T{y) ~ y n 
as y — ► 0, the two- variable scaling function gives a sin- 
gular correction to the power-law near the critical point: 
F(x,r) ~ x~ a F(x c r) - X " a+ " c for x « r" 1 ^: only 
when x ~ r" 1 ^ will the power-law x~ a be observed. 
This is what happened in two dimensions to the inte- 
grated avalanche size distribution (figures |35| and p6| ) and 
the avalanche correlation functions (figure~p7|b). 

For the fourth example, we address finite-size scaling of 
a property K of the system, as we vary a parameter r. If 
we measure K(r, L) for a variety of sizes L (say, all with 
periodic bo und ary conditions), we expect (in complete 
analogy to (|A^)) 



(A10) 



Now, suppose our "thermometer" measuring r is weakly 
size-dependent, so the measured variable is C(r) = r + 
c/L + C2/L 2 + . . . The effects on the scaling function is 



K 



:(c(t),L \ - r~" x 
(cL 1 ^- 1 + c 2 L 1 ' v - 2 ) IC'irL 1 '") + ...). (All) 



In two and three dimensions, v > 1 and these correc- 
tion terms are subdominant. In four and five dimensions, 
we find 1/2 < v < 1, so we sh ould include the term mul- 
tiplied by c in equation (All). However, we believe this 
first term is zero for our problem. For a fixed boundary 
problem (all spins "up" at the boundary) with a first or- 
der transition, there is indeed a term like c/L in r{L) |15| . 
At a critical transition, the leading correction to r(L) can 
be c/L or a higher power of L {l/L 2 and so on). This 
seems to depend on the model studied, the geometry of 
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the system, and the boundary conditions (free, periodic, 
ferromagnetic, . . .) |46| . Furthermore, for the same kind 
of model, the coefficient c itself depends on the geom- 
etry and boundary conditions, and it can even vanish, 
which leaves only higher order corrections. In a periodic 
boundary conditions problem like ours, we expect that 
the correction is smaller than c/L. Our finite-size scaling 
collapses for spanning avalanches N, the second moments 
(S 2 ), and the magnetization jump AM, were successfully 
done by letting c = 0. 

For the fifth example, consider a property K depending 
on two external parameters: r (the disorder for example) 
and h (could be the external magnetic field H — H c ). 
Analogous to (|A9|), K should then scale as 

K(r, h) ~ r'" K(h/r 06 ). (A12) 

Consider now the likely dependence of the field h on the 
disorder r. A typical system will have a measured field 
which depends on the randomness: C{h) = h+br+b2r 2 + 
. . . (Corresponding nonlinearities in the effective value of 
r are subdominant.) This system will have 



F(x, r, h) ~ x^Tix^r, x C!3S h). 



(A15) 



K(r,C(h) 



K-ih/r 136 ) + (br + b 2 r 2 ) r~ 0S 1C'(h/r 05 ) 



(A13) 



Now, for our system 1 < (35 < 2 for dimensions three and 
above. This means that the term multiplied by b is domi- 
nant over the critical scaling singularity: unless one shifts 
the measured h to the appropriate h! = h + br, the curves 
will not collapse (e.g., the peaks will not line up horizon- 
tally). We measure this (non-universal) constant for our 
system (Table §), using the derivative of the magnetiza- 
tion with field dM/dH(r, h). The magnetization M(r, h) 
and the correlation l ength £(r, h) should also collapse ac- 
cording to equation ( A12| ) (but with h + br instead of h) ; 
we don't directly measure the correlation length, and the 
collapse of M(r, h) in figure |l6| b includes the effects of the 
tilt b. In two dimensions, (35 is large (probably infinite), 
so in principle we should need an infinite number of cor- 
rection terms: in practise, we tried lining up the peaks 
in the curves (with no correction terms); because we did 
not know (3 (which we usually obtained from AM, which 
gives (3/v = in two dimensions), we failed to extract 
reliable exponents in two dimensions from dAI/dH. 

For the sixth example, suppose F depends on r, h, and 
a size x. Then from the previous analysis, we expect 



F(x,r,h) ~ x- a F(x C r, h/r^ A ) 



(A14) 



Notice that universality only removes one variable from 
the scaling form. One could in practice do two-variable 
scaling collapses (and we believe someone has probably 
done it), but for our purposes these more general scaling 
forms are used by fixing one of the variables. For exam- 
ple, we measure the avalanche size distribution at various 
values of h (binned in small ranges), at the critical dis- 
order r = 0. We can make sense of equation (A14) by 
changing variables from h/r^ s to x^ s h: 



Before we can set r = 0, we must see what are the pos- 
sible corrections to scaling in this case. If the disorder r 
depends on the field, then instead of the variable r, we 
must use r + ah (the analysis is analogous to the one in 
example five; other corrections are subdominant). Set- 
ting r = now, leaves F dependent on its first variable, 
as well as the second: 



F(x,r,h)~ x~ a F(x<{ah), x^ 5 h) « 
(jT(0, x^ s h) + 

ahxi F^(0, x^ 5 h) 



(A16) 



where F^ 1 ' ^ is the derivative of T with respect to the 
first variable (keeping the second fixed). 

For the binned avalanche size distribution, x^ is S a , 
where 0<cr<l/2as we move from two dimensions to 
five and above. Thus, the correction term will only be 
important for rather large avalanches, S > h~ 1 / a , so long 
as we are close to the critical point. Expressed in terms 
of the scaling variable, important corrections to scaling 
occur if the scaling variable X = S a ^ s h > h}~P 6 . For us, 
(35 > 3/2, and we only use fields near the critical field 
(h < 0.08), so the corrections will become of order one 
when X = 4 for the largest h we use. In 3 and 4 di- 
mensions, this correction does not affect our scaling col- 
lapses, while in 5 dimensions some of the data needs this 
correction. We have tried to avoid this problem (since we 
don't measure our data such that it can be used in a two- 
variable scaling collapse) by concentrating on collapsing 
the regions in our data curves where this correction is 
negligible. 

A similar analysis can be done for the avalanche time 
distribution, which has two "sizes" S and t and one pa- 
rameter r which is set to zero; beca use we integrate over 
the field h the correction in (A16) does not occur, and 
other scaling corrections are small. 

Finally, we discuss the unusual exponential scaling 
forms we developed to collapse our data in two dimen- 
sions. If we assume that the critical disorder R c is zero 
and that the linear term in the rescaling of r vanishes 
(C^er in equation ( A7) vanishes), then from sym metry the 
correction has to be cubic, and equation (A7) becomes: 



F(x,r) = (1-ae) F((l-e)x, (l + /cer 2 )r). (A17) 

with k (which is not universal) and a constants, and e 
small. 

Taking the derivative of both sides with respect to 
e and setting it equal to zero gives a partial differen- 
tial equation for the function F. To solve for F, we 
do a change of variable: (x, r) — > (x, y) with y — 
x e~ a / r . The constant a* is determined by requiring 
that y rescales onto itself to order e: we find a* = 1/2 k. 
We then have: 
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= 



F(x,y) 



dF 
dx 



which gives 



F{x,r) = x- a P{xe- 1 ' 2kr2 y 



(A18) 



(A19) 



This is one of the forms we use in 2 dimensions for 
the scaling collapse of the second moments (S 2 )i n t, the 
avalanche size distribution Di n t integrated over the field 
H, the avalanche correlation Gi n t, and the spanning 
avalanches N. We use another form too which is ob- 
tained by assuming that the critical disorder R c is not 
zero but that the linear term in the rescaling of r still 
vanishes. Instead of equation (A17), we have: 



F(x,r) = (1 - ae) f{(1 - e)x, (l + Zer)rY (A20) 
The function F becomes: 



F(x,r) 



x 



T^xe 



l/lr 



(A21) 



The c orrec tions t o sca ling for the last two forms (equa- 
tions ( A19 ) and flA2l|) ) are similar to the ones discussed 
above. They are all are subdominant. 



APPENDIX B: FULL DERIVATION OF THE 
MEAN FIELD SCALING FORM FOR THE 
INTEGRATED AVALANCHE SIZE 
DISTRIBUTION 

The mean field scaling form for the integrated 
avalanche size distribution Di nt (S, R) was obtained in 
section IV A using the scaling form of the avalanche size 
distribution D(S, R, H). The scaling form for D int (S, R) 
can also be obtained by integrating the avalanche proba- 
bility distribution D(S,t) (derived originally in di- 
rectly: 



The average mean field magnetization M and the 
avalanche probability distribution D(S,t) are given by 

@0: 



M(H, R) = l-2 J 



-JM(H)-H 



p(f) df, 



(B2) 



and 



D(S,t) 



S 



5-2 



(5-1)! 



(t+l) 3 - 1 e- 3 ^ 



(B3) 



To solve equation (|B1[ ), let's define the variable y 
(—JM — _ff)/(v / 2 R) and rewrite the integral as: 



D mt (S,R) = V2R 
[ P (V2Ry) D(S, 2Jp(V2Ry) - l) 

(l - 2 Jp(V2Ry)) dy 

where we have used: 



dy_ = _J_ 

dH ^2R 



2 p(—JM — H) 
~ J 1 - 2Jp(—JM — H) ~ 1 



(B4) 



(B5) 



Since we are interested in the behavior of the integrated 
avalanche distribution for large sizes, the factorial in 
equation (B3) can be expanded using Stirling's formula. 
To first order, we have: 



(5-1)! 



S s V2^ 
e s VS 



(B6) 



Substituting this and the random field distribution func- 
tion p, 



p(V2Ry) = -jL- e- y \ 
f^R 



D, 



V2n 

/+oo 
p(-JM - H) D(S,t) dH (Bl) in equation Q ? we obtain: 



(B7) 



where p(—JM — H) is the probability distribution for 
the random fields, and p(—JM — H) dH is the prob- 
ability for a spin to flip between fields —JM(H) — H 
and -JM(H + dH) - (H + dH). D(S,t) is the prob- 
ability of having an avalanche of size 5, a small "dis- 
tance" t = 2Jp{-JM - H) - 1 from the infinite 
avalanche at p(—JM — H) = 1/2 J, given that a spin 
has flipped at - JM - H ||||]. (The scaling form for 
the non-integrated avalanche size distribution D(S, R, H) 
(eqn.|l9|) is obtained from D(S,t) by expressing t as a 
function of R and H fb|,[l4| ) . J is the coupling of a spin 
to all others in the system, H is the external magnetic 
field, and R is the disorder. The advantage of this proce- 
dure is that we can find out something about the scaling 



+ OC 



Rc 
R 



dy (B8) 



where C = S~i e s R c /(ttRV2), and 5 is large. 

For disorders above but close to the critical disorder 
i? c , we have: 



s - 



function T> 



(int) 



Rc 
R 

Mr) (^Y 



(B9) 
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If we assume that only terms up to 5 (1 — R c /R) 2 are 
important (terms like 5 (1 - R c /R) 3 and 5 (1 - R c /R) 4 
go to zero as i? — > i? c ), and we note that the integrand 
in equation (B8) is an even function of y, equation flB8| ) 
becomes: 



D int (S,R) « 2C x 



1_ if 6 



(BIO) 



The asymptotic behavior of the above integral, as 5 —> 
oo, is obtained using Laplace's method p7j]. The idea is 
as follows. The asymptotic behavior as 5 — > oo of the 
integral: 



I(S) = / f(x) e s ^ x) dx 



(Bll) 



can be found by integrating over a small region [c— e, c+e] 
(instead of the interval [a, 6]) around the maximum of the 
function 4>(x) at x = c, since in the asymptotic expan- 
sion, the largest contribution to the integral will be from 
this region. The corrections will be exponentially small. 
The maximum of cj> must be in the interval [a, b], f(x) 
and 4>{x) are assumed to be real continuous functions, 
and /(c) 7^ 0. f(x) and <p(x) can now both be expanded 
around x = c, and the integral solved. Often the integral 
is easier to handle if the limit of integration is extended 
to infinity. This will add only exponentially small correc- 
tions in the asymptotic limit of S — ► oo. 



Let's apply this method to equation (BIO). The func- 
tion in the exponential has a maximum at y = 0. The 



function [ 1 — ^ 



is not zero there if R ^ R c . We 



can thus expand both functions in the integral of equa- 
tion (BIO) around y — 0. Defining u = y 2 V~S, we obtain: 



D int (S,R) « C x x 



~ ~R J + ~R 7f ~ 2~R S 7 " 



\fu 



where 



Ci = 



irV2 R 



R c Q _| - 
A 4 e 



(B12) 



(B13) 



5 is large, i? is close to but not equal to i? c , and only 
terms up to 5(1 — R c /R) 2 are non- vanishing. In the 
asymptotic limit of 5 — > oo we can ignore terms with 
powers of 5 in the denominator, and look at the distri- 
bution for R close to R c . To first order in r = (R c —R) / R, 
R c /R rs 1 and 1 — R c /R « — r, which gives: 



D int (S, R) 



e\ J (-rVS + u) — (B14) 



where we have expanded the integration to infinity. As 
mentioned above, this will only add exponentially small 
corre ction s in the asymptotic limit of 5 — > oo. Equa- 
tion (B14) is the integrated avalanche size distribution 
in mean field for large sizes 5, and finite Sr 2 . We see 
right away that it gives the correct scaling form: 



D int (S, R) 



V 



(int 



(Vs 



(B15) 



where ± indicates the sign of r, the exponent T + af35 and 
a are 9/4 and 1/2 respectively, and the scaling function 



V { l nt) is: 



V 



(int 



b/JL M) 2 



f±(Vs 



(B16) 



The function T±(\/~S \r\j is proportional to the integral 
in equation (B14). Note that the above result is equiv- 



alent to the one obtained (eqn. |23|) by integrating the 
scaling form of D(S, R, H) over the field H . 

What is the behavior of the scaling function X^™*' (X) 

for small and large positive ar gume nts X — >/ S (— r) > 
(R> R c )l From equations pOj ) and ( p!6| ), for small 
arguments we have a polynomial in X: 



V { l nt) {X) 



A + BX + CX 2 + 0(X 3 



(B17) 



These parameters can be calculated numerically. We ob- 
tain in mean field: 



V 



(int) 



0.232 + 0.243X-0.174X 2 
0.101X 3 + 0.051X 4 



On the other hand, for large arguments we find: 



^ nt) (X)«7r 1/2 e-^ Vx(l + 0(X- 



(B18) 



(B19) 



In general, for all dimensions, in equation ( p319| ) the ex- 
ponential is of X 1 /® (1/cr = 2 in mean field), since the 
exponent a gives the exponential cutoff to the power law 
distribution for large X, and the power of X is /3 (/3 = 1/2 
in mean field). One can see the latter by expanding the 
distribution function Di nt (S,R) in terms of 1/5 (5 is 
large), analogous to pgfl: 



D int (S,R) = J2f n (r) S~ 



(B20) 



n=l 



Since X = S a (~r), then we can write 5 = X 1 ^ (-r)- 1 / 
and obtain: 



3G 



D int (S,R) = J2f n (r) X- n '°{-ry 



(B21) 



71=1 



The scaling function V K _>{X) scales like S*( r+(J ' 3<5 ) 
D int {S, R): 

V { " lt) {X) r 



Y,fn{r) x- n / a X ( - T +' yf3S '>/' J x 

(_ r )«/°'(_ r )-(- r + cr /3< 5 )/ CT 



n=l 



(B22) 



and since it is only a function of X , it must satisfy: 



T> { * nt) (X)~ J2 5" X- n ^X^ T+a ^/' 7 (B23) 

n=l 

where g n is independent of r. 

The exponent combination (r + uf35) j a can be rewrit- 
ten as: 

I±£gg = 2 + r + ^-2 = 2 

(T(T(T (7 

where we have used the scaling relation |l^Jl^]: /3 — 
(38 = [t — 2)/cr. Thus we have for the scaling function 

V { l nt) {X): 



■D { " lt) (X) ~ X -9« X ~ n/a = xP £(X 1/a ) (B25) 
n=-l 



which shows (compare to equation ( B1S| )) that the power 
of X is indeed the exponent (3. 

We have used the results of the expansion of the mean 
field scaling function f>^ n t '{X) f or sm all and large pa- 
rameters (equations (B17) and ( B 1 9| ) ) , to build a fit- 
ting function to the integrated avalanche size distribution 
scaling functions in 2, 3, 4, and 5 dimensions, described 
in section IV B. 



Finally, note from equation (B17) that the scaling func- 
tion: 



S m f is the total size of the system. We want to derive 
the scaling form for the number of such avalanches in 
half of the hysteresis loop (for H from — oo to +oo) as 
a function of the system size S m f and the disorder R. 
The number of mean field spanning avalanches is pro- 
portional to the probability of having avalanches of size 
larger than y/S m f. Since we want the number of span- 
ning avalanches, we need to multiply this probability by 
the total number of avalanches. For large system sizes, 
this scales with the system size S m f (corrections are sub- 
dominant). We thus obtain by integrating over equation 
( pT5| ) (which gives the scaling form for the probability 
distribution of avalanches of size 5* in the hysteresis loop): 




S~* e 



N m f(S m f, R) 

</S|r|) 



•F±(VS|r|) dS 



(CI) 



Let's define u — >/S\r\, then equation (CI) can be writ- 
ten as: 



N mf (S mf ,R) ~ 2 S mf |r|* x 



i u 2 e 2 T±(u) du 



(C2) 



The integral X is a function of S , ^ i y|'"| only, and we can 
write it as: 



l=\S mf \r 



1 N? f (S, 



7 mf\ 



(C3) 



to obtain the scaling form for the number N m f of mean 
field spanning avalanches: 

N mf (S mf ,R) ~ Sl f M? f (sl f \r\) (C4) 

From this scaling form, we can extract the exponents 
6 = 3/8, and 1/v = 1/4. This form is used for collapses 
of the spanning avalanche curves in mean field (see mean 
field section). 



(B26) 



used earlier in reference |l3j], is not analytic for small 
arguments Sr 2 , from which we conclude that the appro- 
priate scaling variable should be VS* (— r) and not Sr 2 . 
(Notice that this no longer seems true in two dimensions; 
sec section on 2 dimensional results.) 



APPENDIX C: DERIVATION OF THE MEAN 
FIELD SCALING FORM FOR THE SPANNING 
AVALANCHES 

We have defined earlier a mean field spanning 
avalanche to be an avalanche larger than yj S m f, where 
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3d 


4d 


5d 


mean field 


R c 


2.16 ±0.03 


4.10 ±0.02 


5.96 ±0.02 


0.79788456 


H c 


1.435 ±0.004 


1.265 ±0.007 


1.175 ±0.004 





b 


0.39 ±0.08 


0.46 ±0.05 


0.23 ±0.08 






TABLE I. Numerical values for the critical disorders and 
fields, and the "tilt" parameter b (see section on magnetiza- 
tion curves collapses) in 3, 4, and 5 dimensions extracted from 
scaling collapses. The critical disorder is obtained from col- 
lapses of the spanning avalanches and the second moments of 
the avalanche size distribution. The critical field is obtained 
from the binned avalanche size distribution and the magne- 
tization curves. H c is affected by finite sizes, and systematic 
errors could be larger than the ones listed here. The mean 
field values are calculated analytically Jl|[14|. The "tilt" b 
is obtained from the dM/dH collapses using the values for 
the critical disorder and field from this table and the values 
for the exponents from Table 



III, Only the parameter b 



allowed to vary. The values in 2 dimensions (which are not 
listed) are less accurate. Depending on the scaling form we 
obtain a critical disorder of 0, 0.45, or 0.54. The critical field is 
around 1.32 and is estimated from the binned in H avalanche 
size distribution and magnetization curves (see text). The 
"tilt" b was not measured. R c , H c , and b are not universal 
characteristics of the system. 



39 



measured exponents 3d 4d 5d mean field 

TJu 0.71 ±0.09 1.12 ±0.11 1.47 ±0.15 2 

6 0.015 ±0.015 0.32 ±0.06 1.03 ±0.10 1 

(r + a/38 -3)/av -2.90 ±0.16 -3.20 ±0.24 -2.95 ±0.13 -3 

l/a 4.2 ±0.3 3.20 ±0.25 2.35 ±0.25 2 

r + a/38 2.03 ±0.03 2.07 ±0.03 2.15 ±0.04 9/4 

t 1.60±0.06 1.53 ±0.08 1.48 ±0.10 3/2 

d + j3/v 3.07 ±0.30 4.15 ±0.20 5.1 ±0.4 7 (at d c = 6) 

0/u 0.025 ±0.020 0.19 ±0.05 0.37 ±0.08 1 

avz 0.57 ±0.03 0.56 ±0.03 0.545 ±0.025 1/2 



TABLE II. Values for the exponents extracted from scaling collapses in 3, 4, and 5 dimensions. The mean field values 
are calculated analytically p^ , ^4[ . v is the correlation length exponent and is found from collapses of avalanche correlations, 
number of spanning avalanches, and moments of the avalanche size distribution data. The exponent 6 is a measure of the 
number of spanning avalanches and is obtained from collapses of that data, (r ± a (38 — 3)/av is obtained from the second 
moments of the avalanche size distribution collapses. 1/a is associated with the cutoff in the power law distribution of avalanche 
sizes integrated over the field H, while r ± a/38 gives the slope of that distribution, r is obtained from the binned avalanche 
size distribution collapses. d + /3/v is obtained from avalanche correlation collapses and f3/v from magnetization discontinuity 
collapses, avz is the exponent combination for the time distribution of avalanche sizes and is extracted from that data. 
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calculated exponents 3d 4d 5d mean field 

a05 0.43 ±0.07 0.54 ±0.08 0.67±0.11 3/4 

05 1.81 ±0.32 1.73 ±0.29 1.57 ±0.31 3/2 

0.035 ±0.028 0.169 ±0.048 0.252 ±0.060 1/2 

av 0.34 ±0.05 0.28 ±0.04 0.29 ±0.04 1/4 

V = 2 ± (0 - 05) /v 0.73 ±0.28 0.25 ±0.38 0.06 ±0.51 



TABLE III. Values for exponents in 3, 4, and 5 dimensions that are not extracted directly from scaling collapses, but instead 
are derived from Table |n]and the exponent relations (see [ p^Jlrj ]). The mean field values are obtained analytically Both 
a05 and 05 could have larger systematic errors than the errors listed here. See the binned avalanche size distribution section 
for details. 
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Exponents 




L=10,20,40 




L=20,40,80 


and R c 


r = (R c -R)/R 


r = (Rc - R)/R c 


r = (R c - R)/R 


r = (R c - R)/R c 


l/v 


0.96 ±0.07 


1.07 ±0.05 


1.05 ±0.10 


1.12 ±0.06 





0.35 ±0.10 


0.34 ±0.06 


0.32 ±0.04 


0.32 ±0.04 


R c 


4.09 ±0.02 


4.09 ±0.01 


4.095 ±0.015 


4.10±0.01 



TABLE IV. Exponent values and critical disorder R c from collapses of spanning avalanche curves in 4 dimensions. Three 
curves (different linear size L) are collapsed together, with r = (R c — R)/R and r — (R c — R)/R c . Tables [iv|, and ^ give 
information equivalent to that given in for example figures []a and Graphs showing two points with an extrapolation to 
1/L seemed unnecessary. 
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Exponents 


L=10,20,40 




L=20,40,80 




r = {R c -R)/R r = (R c -R)/R c 


r={R c - R)/R 


r = (R c - R)/R c 


l/v 


1.10 ±0.04 1.24 ±0.08 


1.10 ±0.05 


1.11 ±0.05 


P/u 


0.195 ±0.035 0.19 ±0.05 


0.18 ±0.05 


0.20 ±0.06 



TABLE V. Exponent values for \ jv and P/v, obtained from scaling collapses of the change of the magnetization AM due 
to the spanning avalanches. Three curves of different size L are collapsed together with r = (R c — R)/R and r — (R c — R)/R c , 
where R c = 4.10 ± 0.02. See also comment in Table [F^. 
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Exponents 


L= 


=5,10,20 


L= 


10,20,30 




r = (R e -R)/R 


r=(R c - R)/R c 


r = (R c - R)/R 


r = (Rc - R)/Rc 


1/v 


1.40 ±0.05 


1.60 ±0.10 


1.41 ±0.07 


1.53 ±0.13 


-(r + af3S-3)/av 


2.75 ±0.10 


2.70 ±0.10 


2.93 ±0.08 


2.90 ±0.08 



TABLE VI. Exponent values from the collapses of second moments of the avalanche size distribution curves in 5 dimensions. 
Three curves of different size L are collapsed together with r = (R c — R) / R and r — (R c — R) / R c , where R c — 5.96 ±0.02. See 



also comment in Table IV 
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2d 


1/ V 


(r + afiS - 3)/av 


a 


t + a/38* 


T 


av 




conj. 
meas. 




0.13 ±0.13 


-2 

-1.9±0.1 




0.10 ±0.02 


2 

2.04 ±0.04 


3/2 
0.0 ±0.0 


1/2 
0.51 ±0.08 




0.03 ±0.06 



TABLE VII. Conjectured and measured values for some exponents in 2 dimensions. We don't have a conjectured value for 
the exponent combination ouz, but the measured value is 0.64 ±0.02. (*) Note that the distribution of avalanche sizes at R c in 
two dimensions, integrated over the loop, will scale as g-( T +< T l 3S )+'^ ] where the correction uj ~ 1 is due to the singularity in the 
scaling function V i f nt \X) ~ X w as X — > 0. See text section IV, figure^, and appendix A for details. A similar argument can 
be made for the avalanche correlation measurement (integrated over the field H), where due to the singularity of the scaling 
function Q-, the scaling for small x\r\" is x ~ ( - d+l3 ^' l+u ' , with Co ~ 1 (see text and figure ^b). 



2 dimensions 


e 


af3S 


06/v 


av 


1/8 


V V 


conjectured 





1/2 


1 


1/2 





1 2 



TABLE VIII. Conjectured values for some exponents in 2 
dimensions. These exponents were not extracted from col- 
lapses (see text). 
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